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ABSTRACT 


Consider  a given  homogeneous  or  inhomogeneous  linear  difference 

equation  Z.g-o  “ g(r)  where  i ^2  and  r = 0,1,2,...  . 

Suppose  y is  a solution  of  this  equation  and  u,v  are  solutions 
of  the  homogeneous  form  of  this  equation  such  that  u(r)/v(r)  ->■  0, 
y(r)/v(r)  -*•  0,  u(r)/y(r)  0.  It  is  knox^n  that  under  these  circum- 

stances algorithms  for  the  computation  of  y based  on  forward 
recurrence  or  backward  recurrence,  such  as  the  Miller  algorithm, 
are  numerically  unstable". 

Stable  algorithms,  such  as  the  method  of  Olver  in  the  case  1=2, 
have  been  based  on  approximating  y(r)  by  the  solutions  of  a 

certain  sequence  of  boundary  value  problems.  More  specifically,  Y^^^) 

is  a solution  of  the  difference  equation  that  coincides  with  y(r) , 
over  some  fixed  initial  range  of  r,  say  r = i,  i+1, . . . , i+j -1 , and 
satisfies  Y^Cr)  = 0 for  r = n,n+l, , . . ,n+£-j-l.  Here  j is  an 

integer  whose  value  depends  on  the  asymptotic  behavior  of  the  chosen 
solution  y(r)  and  n is  an  arbitrary  large  integer.  Boundary  value 
problems  of  this  type  are  shown  to  be  equivalent  to  tx<ro  initial  value 
problems  of  order  j and  £-j  by  factorization  of  the  linear  differ- 
ence operator.  The  solution  of  the  problem  of  order  j is  obtained 
by  forward  recurrence;  the  solution  of  the  other  problem  is  obtained 
by  backward  recurrence. 

The  algorithm  is  specified  completely  for  a broad  class  of  linear 

difference  operators.  This  class  includes,  for  example,  every  constant 

coefficient  operator.  Convergence  of  v (r)  to  yCr)  as  n for 

^ n 

fixed  r is  proved  and  an  expansion  of  the  truncation  error  is  derived 
Numerical  stability  is  demonstrated  under  appropriate  conditions.  The 
method  is  tested  by  numerical  examples  involving  fourth-order  equations 
with  variable  coefficients. 
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CHAPTER  1.  INTRODUCTION 


1 . 1 Miller's  and  Olveras  Algorithms 

Linear  recurrence  relations  (linear  difference  equations)  are  satis- 
fied by  many  of  the  higher  transcendental  functions  of  mathematics,  and 
they  arise  quite  naturally  in  applications  of  mathematics.  Discretization 
of  ordinary  differential  equations  in  numerical  analysis  is  one  important 
area  in  which  they  occur.  Because  of  their  special  form,  linear  recurrence 
relations  are  especially  well  suited  to  the  construction  of  computing 
algorithms.  Stability  requirements  in  the  case  of  second-order  homogeneous 
relations  led  to  the  well-known  Miller  backward  recurrence  algorithm  [ 5 ] , 
which  was  developed  for  the  computation  of  modified  Bessel  functions.  This 
reference  was  published  in  1952.  The  review  paper  of  Gautschi  [9]  pro- 
vides a survey  of  algorithms  and  applications  in  the  case  of  second-order 
homogeneous  equations.  Nowadays,  the  name  of  Miller  often  attends  when  a 
minimal  solution  (defined  below)  of  a homogeneous  equation  of  arbitrary 
order  is  computed  using  backward  recurrence  [28].  These  algorithms  always 
involve  backward  recurrence  from  assumed  starting  values  (since  these  are 
not  known  in  general)  followed  by  a normalization  procedure,  such  as  match- 
ing one  value  of  the  computed  solution  with  the  desired  solution  in  order 
to  arrive  at  the  proper  scale  factor  which  relates  these  two  solutions. 

A different  class  of  algorithms  becomes  necessary  when  "intermediate" 
solutions  of  a linear  recurrence  relation  are  to  be  computed.  Intuitively 
speaking,  an  intermediate  solution  is  one  for  which  two  other  solutions  can 
be  found  such  that  one  dominates,  and  the  other  is  dominated  by,  the  inter- 
mediate solution.  We  are  speaking  here  of  dominance  in  the  following  sense: 


if  a = (a(0) , a(l) , . . . ) and  b = (b(0) ,b(l) , . . .)  are  two  infinite 
sequences  then  b dominates  a whenever  a(r)/b(r)  ->•  0 as  r . 

Neither  forward  recurrence  nor  the  Miller  algorithm  is  numerically 
stable  for  an  intermediate  solution.  This  is  because  recurrence  in  either 
direction,  when  carried  out  in  finite  precision,  will  be  contaminated  by 
small  components  of  the  other  solutions.  An  analysis  of  the  prop- 
agation of  these  rounding  errors  appears  in  [18]  . In  the  forward  direction, 
as  the  recurrence  proceeds  through  more  and  more  steps,  the  dominant  solu- 
tions grow  and  eventually  overwhelm  the  desired  intermediate  solution. 

Since  our  concept  of  dominance  is  really  dominance  "at  infinity",  the  situa- 
tion for  recurrence  in  the  backward  direction  is  not  completely  analogous. 
However,  if  we  start  at  a sufficiently  large  value  of  r and  recur  backward, 
other  solutions  grow  at  a faster  rate  than  the  wanted  intermediate  solution 
as  r decreases. 

In  addition  to  intermediate  solutions,  a difference  equation  has  minimal 
(or  "recessive")  and  maximal  (or  "dominant")  solutions.  Intuitively,  a 
minimal  solution  does  not  dominate,  and  a maximal  solution  is  not  dominated 
by,  any  other  solution  of  the  difference  equation  as  r . Ultimately, 

forward  recurrence  is  stable  for  maximal  solutions  since  the  error  cannot 
grow  more  rapidly  than  the  solution.  Similarly,  backward  recurrence  is 
stable  for  minimal  solutions,  at  least  in  a region  of  sufficiently  large  r . 


We  shall  regard  the  solutions  of  a linear  difference  equation,  as  well  as 
the  coefficients  of  the  equation,  as  complex-valued  functions  of  a discrete 
real  variable.  We  shall  assume  that  this  discrete  real  variable  proceeds 
in  unit  increments;  thus  the  solutions,  as  well  as  the  coefficients,  of 
linear  difference  equations  are  complex  sequences.  We  shall  be  interested 
only  in  those  linear  difference  equations  whose  coefficient  sequences  are 
of  infinite  length. 
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The  lowest  order  of  difference  equation  for  which  an  intermediate  solu- 
tion can  exist,  and  then  only  when  the  equation  is  inhomogeneous,  is  two. 
Such  an  equation  is  of  the  form 

(1.1.1)  a(r)y(r-l)-b(r)y(r)+c(r)y(r+l)  = d(r)  , r = 1,2,...  . 

Since  the  general  solution  is 

y(r)  = au(r)  + pv(r)  + h(r)  , r = 0,1,2,...  , 

where  a and  p are  constants,  u=  ( u(0) , u(l) , . . . ) and  v=  (v(0),v(l), 
...)  are  linearly  independent  solutions  of  the  homogeneous  equation,  and 
h = (h(0) ,h(l) , . . . ) is  a particular  solution  of  (1.1.1),  there  can  be  at 
most  three  asymptotically  distinct  types  of  behavior  of  solutions  as 
r . Let  us  assume  that  v dominates  u and  that  h is  an  inter- 

mediate solution.  In  addition,  assume  u(0)  ^ 0 and  v(r)  0 for  all 
sufficiently  large  r . Under  these  assumptions  a stable  algorithm  for 
computing  intermediate  solutions  of  (1.1.1)  was  presented  by  Olver  in  [21]. 
If  y = (y (0) , y (1) , . . . ) is  the  desired  intermediate  solution.  Giver's 

algorithm  approximates  y by  a sequence  of  solutions  y = (y  (0),y  (1), 

n n n 

...)  such  that  y (0)  = y(0)  and  y (n)  = 0 where  n is  an  appropriate- 

n n 

ly  chosen  positive  integer.  Note  that  only  one  initial  value  of  the  wanted 
intermediate  solution  is  required  to  determine  the  sequence  of  approxima- 
tions. The  sequence  is  shown  to  converge  pointwise  to  y as  n . 

Giver's  algorithm  efficiently  produces  the  first  n-1  unknown  terms 

of  y as  the  solution  of  a tridiagonal  linear  system  of  order  n-1  . The 
n 

procedure  used  for  solving  these  linear  equations  is  equivalent  to  Gaussian 
elimination  without  pivoting.  It  consists,  in  effect,  of  a forward  elimina- 
tion stage  followed  by  a back  substitution  stage.  The  algorithm  is  construct- 
ed so  that  the  numerical  stability  of  this  procedure  may  be  assessed.  The 


fo-n^ard  elimination  stage  is  realized  by  two  forward  recurrences.  One  of 
these  uses  the  homogeneous  form  of  (1.1.1)  to  generate  a dominant  solution 
P ~ (p(0) >p(l) 5 • • •)  . The  initial  values  used  for  starting  this  recurrence 
are  p(0)  = 0 and  p(l)  = 1 , and  the  solution  p is  asymptotically  pro- 
portional to  V . Therefore,  the  rounding  errors  ultimately  propagate  in 
proportion  to  p , although  they  could  propagate  faster  than  p at  first 
(if  p is  essentially  proportional  to  u for  early  values  of  r , 
corresponding  to  an  excessively  small  value  of  u(0)).  A combined  algorithm 
that  employs  both  Miller's  algorithm  and  Giver's  algorithm  in  adjacent 
ranges  of  r is  proposed  in  [22]. 

The  second  forward  recurrence  in  the  forward  elimination  stage  of 
Giver's  algorithm  uses  the  inhomogeneous  first-order  equation 

a(r)e(r-l)  - c(r)e(r)  = d(r)p(r)  , r = 1,2,..., 

where  the  solution  e = (e (0) , e(l) , . . . ) satisfies  e(0)  = y(0)  . Here 
two  distinct  asymptotic  forms  of  behavior  may  exist,  but  there  is  no  reason 
to  expect  that  the  calculation  of  e will  be  asymptotically  unstable  in 
general.  However,  as  was  the  case  for  the  solution  p , there  may 
be  applications  where  rounding  errors  grow  more  rapidly  then  e(r)  over 
some  finite  initial  range  of  values  of  r . 

The  back  substitution  consists  of  backward  recurrence  of  the  inhomo- 
geneous first-order  equation 

p(r+l)y  (r)  - p(r)y  (r+1)  = e(r)  , r = n-1 , n-2 , . . . , 1 , 
n n 

starting  with  ~ ^ stability  of  this  backward  recurrence  is 

guaranteed  in  the  "asymptotic  region"  belonging  to  the  chosen  solution 
y(r)  , that  is,  the  infinite  range  of  r which  is  such  that  the  asymptotic 
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behavior  of  y(r)  relative  to  u(r)  and  v(r)  is  maintained. 

Giver's  method  of  solving  the  tridiagonal  system  leads  to  a very  con- 
venient representation  of  the  truncation  error.  For  each  r and  for  all 
n > r , the  expansion 


(1.1.2) 


y(r)  = y (r)  + p(r) 
n 


y e(s) 

^ p(s)p(s+i) 
s=n 


is  valid.  Consequently,  for  each  r and  n > r and  v > 1 we  have 


n+v 


n+v-1 

(r)  - y^(r)  = p(r)  J 


e(s) 


s=n 


p(s)p(s+l) 


This  expansion  is  useful  as  a basis  for  deciding  on  an  appropriate  value  of 
n , i ■ e . , for  deciding  when  to  stop  the  forward  elimination  and  begin  the 
back  substitution.  Suppose  that  for  particular  values  of  e > 0 and 
r , say  r = m , it  is  desired  that 

ly(m)  - y^(m) I < s . 


A value  of  v(>  1)  is  selected  and  yfm)  is  approximated  bv  y , (m)  . 

■ n+v 

The  problem  then  is  to  determine  the  smallest  value  of  n such  that 


< £ . 


This  is  done  by  continuing  the  forward  elimination  until  r = m at 
which  point  calculation  of  the  numbers 


r\(r) 


r+v-1  . 

!p(ni)  I , T N I . r = m,m+l,... 

p(s)p(s+l)' 


is  begun.  Forward  elimination  is  continued  further,  together  with  the 
calculation  of  'n(^)  > until  a value  of  r is  reached  such  that  Ti(r)  < e . 
This  value  of  r is  taken  for  n ; its  adequacy  may  be  checked  by  taking 
additional  terms  in  the  infinite  expansion  of  y(m)  - Y^(^)  » given  by 

(1.1.2)  with  r = m . 

In  addition  to  being  appropriate  for  computing  intermediate  solutions 
of  (1.1.1),  Giver’s  algorithm  is  suited  to  computing  the  minimal  solution 
of  a homogeneous  equation.  It  can  therefore  be  compared  to  Miller's 
algorithm.  The  advantage  of  Giver's  algorithm  over  Miller's  algorithm 
is  that  the  latter  does  not  have  a criterion  for  determining  automatically 
an  appropriate  value  for  n . Gn  the  other  hand,  if  an  adequate  value 
of  n is  known  from  other  considerations  (as,  for  example,  from  an 
asymptotic  estimate),  then  Miller's  algorithm  may  require  less  computa- 
tional effort. 

Zahar  [27],  Gliver  [19,20],  and  Cash  [ 6]  have  considered  extensions  of 
Giver's  algorithm  to  linear  difference  equations  of  higher  order.  In  the 
present  work  we  shall  construct  and  analyze  an  algorithm  which  has  features 
in  common  with  these  extensions.  We  shall  also  illustrate  the  algorithm 
by  means  of  numerical  examples. 
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1.2  Factorization  of  Linear  Difference  Equations 

Let  F be  the  field  of  complex  numbers  and  let  S be  the  set  of  all 
complex-valued  infinite  sequences  x = (x(0) ,x(l) , . . . ) . Clearly  S is 
an  infinite-dimensional  linear  space  over  F . Let  be  the  set  of 

all  infinite  upper- triangular  band  matrices  of  the  form 

dp(0)  ....  d^(0)  0 

d„(l)  ....  d^(l) 

dg(r)  ....  d^(r) 

0 

where  d^,d^,...,d^  € S and  d^Cr)  ^ 0 , d^(r)  0 for  infinitely 

many  values  of  r . We  shall  also  use  the  concise  notation  [d^ j d^ , . . . , d^] 
to  represent  infinite  matrices  in  the  set  . If  we  regard  the  infinite 

sequences  in  5 as  infinite  column  vectors,  then  an  infinite  sequence 
Dx  € S is  defined  for  each  D ^ , x € S by  ordinary  matrix  multiplica- 

tion. In  fact,  the  terms  of  Dx  are  given  by 

£ 

Dx(r)  = y d (r)x(r+s)  , r = 0,1,2,...  . 
s=0  ® 

Clearly,  each  D € is  a linear  operator  on  S under  this  definition. 

We  shall  say  that  is  the  class  of  linear  difference  operators  of  order 

£ on  S . 

Note  that  is  defined  for  each  integer  £ > 0 . We  pointed  out  in 

§1.1  that  intermediate  solutions  of  linear  difference  equations 


can  exist  only  if  the  order  of  the  equation  is  at  least  two.  But  we  are 
going  to  factor  linear  difference  operators  into  a product  of  two  lower- 
order  operators.  Therefore,  operators  of  order  one  will  arise.  Further- 
more, the  formalism  to  be  introduced  will  allow  even  the  degenerate  case 
of  operators  of  order  zero. 

If  d^  and  d^  have  no  terms  equal  to  zero,  we  shall  say  the  differ- 
ence operator  [d^ , d^ , . . . , d^]  € is  nonsingular;  compare  [14,  §12.0]. 

The  terms  of  d^  and  d^  will  be  called  the  leading  and  trailing  coeffic- 
ients of  [d^,d^, . . . ,d^]  , respectively. 

We  shall  use  the  notation  to  indicate  a sequence  from  S whose 

first  term  is  x(i)  . Thus 

x^  = (x( i) , x( i+1 ),...)  . 

The  point  i will  be  called  the  initial  point  of  the  sequence.  When  we  do 
not  indicate  a superscript  it  is  to  be  understood  that  zero  is  the  initial 
point  of  the  sequence.  Similarly,  we  shall  use  the  notation  x^’^  to 
indicate  the  (finite)  subsequence  of  x^  whose  final  term  is  x(n)  . Thus 

x^’^  = (x(i) ,x(i+l) , . . . ,x(n) ) . 

The  point  n will  be  called  the  terminal  point  of  the  subsequence. 

Similarly,  for  linear  difference  operators  from  , we  shall  use  the 

notation  or  [d^jdj^,  , . . ,d^]  to  indicate  that  i is  the  initial  point 

of  each  of  the  £ + 1 sequences  involved  in  the  definition  of  the  opera- 
tor. See  Figure  1.  When  the  superscript  is  omitted,  it  is  understood 
that  these  initial  points  are  all  zero  unless  it  is  clear  from  the  context 
(as  in  equation  (1.2.1)  below)  what  these  initial  points  are. 
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Figure  1. 


The  operator  D 


Figure  2.  The  matrix 


Figure  3. 


The  macrix 
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Two  further  notational  conventions  to  be  used  in  this  thesis  are 
Dj  , the  infinite  matrix  obtained  from  by  deleting  the  first  j 

columns  (05j5-£.)  ; and  leading  principal  submatrix  of  order 

n - i + 1 of  Figures  2 and  3.  Note  that  each  diagonal  of 

below  the  principal  diagonal  is  a subsequence  of  one  of  the  sequences 
d^jd^, . . . ,d^_^  j the  sequence  of  d^  is  on  the  principal  diagonal,  and 
the  sequences  d^_^^,...,d^  form  the  diagonals  above  the  principal  dia- 
gonal. The  form  of  the  finite  matrix  D^’'^  is  similar,  with  the  finite 
subsequence  d^’^  lying  on  the  principal  diagonal. 

A linear  difference  equation  is  an  equation  of  the  form  Dx  = g , 
where  g ^ S and  D ^ are  prescribed.  Either  x C S (or  a sub- 
sequence of  X of  the  form  x^’'^)  is  to  be  determined.  This  equation 
is  equivalent  to  the  infinite  linear  algebraic  system 


Dx(r)  = g(r)  , r = 0, 1, . . . . 

Analogously  to  linear  differential  equations,  auxiliary  conditions  are 
required  in  order  to  specify  a unique  solution.  In  our  investigations 

these  auxiliary  conditions  will  consist  of  tv/o  finite  subsequences 

i,i+j-l  , n,n+k-l  . ^ , . , . , 

X and  X (or  possibly  only  one  of  these;  which  are  to  be 

specified  in  advance,  where  j + k = ^ and  j > 0 , k > 0 . These  are 

called  the  initial  and  terminal  conditions , respectively,  when  n > i + j 


mine  y 

(1.2.1) 


Let  us  consider  the  following  finite  boundary  value  problem: 
i+j ,n-l 


Deter- 


such  that 

Dy(r)  = g(r) 


r = i, i+1, . . . ,n-j-l 


together  with  initial  conditions 


v(i+r)  = a , r = 0,l,...,j-l 
r 


(1.2.2) 


and  terminal  conditions 


(1.2.3)  y(n+r)  ~ ^ ~ 0,1,..., k-1 

are  satisfied,  where  DcP^,'2.>l,g^-S,jfc  {0,1,...,-^}  , k - > 

i>0,n>i+j  and  ^ F . In  view  of 

(1.2.2)  and  (1.2.3),  the  finite  subsequences  ^ and  ^ 

are  known.  Of  course,  if  j = 0 then  ^ is  null.  Similarly  if 

k = 0.  These  cases  correspond  to  the  absence  of  either  (1.2.2)  or  (1.2.3). 
The  unknoxvn  subsequence  ^ satisfies  a linear  system  of  algebraic 

equations  which  can  be  expressed  in  matrix  terms  using  the  notation  we 
introduced  above.  The  system  is 

(1.2.4) 

i,n-j-l  i+j,n-l  ^ i,n-j-l 
j / S 

provided  that  n - i - j > max(j,k)  . The  first  partitioned  column  vector 
on  the  right  side  of  (1.2.4)  has  j entries  in  the  upper  part,  and  the 
second  partitioned  column  vector  has  k entries  in  the  lower  part.  The 
length  of  the  other  column  vectors  appearing  in  (1.2.4)  is  n - i - j . 

This  is  what  gives  rise  to  the  restriction  n - i - j > max(j,k)  . This 
restriction  is  irrelevant  for  the  applications  we  have  in  mind,  since  we 
shall  be  considering  sequences  of  finite  boundary  value  problems  as  n ->  «>  . 
Note  chat  the  number  of  zeros  represented  by  the  symbol  0 in  the  parti- 
tioned column  vectors  is  n - i - 2j  in  the  first  vector  and  n - i - 
in  the  second  vector.  Obviously,  the  partitioned  column  vectors  reflect 
the  influence  of  the  boundary  conditions  (1.2.2)  and  (1.2.3). 

Solution  of  the  boundary  value  problem  (1.2.1)  to  (1.2.3)  is  equiva- 
lent CO  solution  of  the  linear  system  (1.2.4).  Equation  (1.2.4)  possesses 


0 ^ 

0 

0 

^n--d,n-i-l  n,n+k-l 

y 
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a unique  solution  if,  and  only  if,  the  matrix  ^ ^ is  nonsingular. 

If  in  addition  we  assume  that  the  difference  operator  D is  nonsingular, 
then  the  solution  ^ of  (1.2.1)  to  (1.2.3)  extends  uniquely  to 

y € S by  forward  and  backward  recurrence.  However,  nonsingularity  of  D 
is  not  necessary  for  this  unique  extension  to  exist;  for  example,  some 
of  the  leading  and  trailing  coefficients  of  D which  appear  in  the 

• • * -r.  ^ 5 n™"  T " 1 . , , 

finite  matrix  D.  could  be  zero. 

3 

Let  us  consider  the  nonsingular  triangular  cases  of  (1.2.4).  The 
first  is  given  by  j = -^  . Then  (1.2.4)  assumes  the  form 


(1.2.4a) 


•ci-t-l  i+£,n-l 

y 


i, n-£-l 
g 


^i,i+^-l  i,i+Z-l 
°0  ^ 


Writing  out  these  equations  in  detail,  we  have 

1-1 

d.(i)y(i+£)  = g(i)  - \ d (i)y(i+s) 

^ s=0  ^ 


1-1 

d^_^(i+l)y(i+^)  + d^(i+l)y (i+^+1)  = g(i+l)  - \ d^ (i+l)y (i+s+1) 

s=0 


I 

\ d (i+£)y(i+^+s)  = g(i+^) 
s=0  ® 


\ d (n--^-l)y(n-^-l+s)  = g(n-£-l) 
s=0  ® 

Since  this  system  is  lower  triangular,  it  is  solvable  by  determining 
y{±+V)  from  the  first  equation,  y(i+£+l)  from  the  second  equation,  and 
so  on.  This  process  is  equivalent  to  developing  y (i+^) , y (i+£+l) , . . . , y (n-1) 
from  (1.2.1)  by  forward  recurrence  starting  with  the  given  values  (1.2.2). 


Then  (1.2.4)  assumes 


The  second  triangular  case  is  given  by  j = 0 . 
the  form 


i,n-l  i,n-l  i,n-l 

•o  y = g 


Writing  out  these  equations  in  detail,  we  have 

i 

I d (i)y(i+s)  = g(i) 
s=0  ® 


n-£,n-l  n,n+^-l 

y 


I 

^ d ^ (n-£-l)y  (n--£-l+s)  = g(n-^-l) 
s=0  ^ 


I 

dQ(n-2)y (n-2)  -f-  d^(n-2)y (n-1)  = g(n-2)  - I d^ (n-2)y (n-2+s) 

s=2 

£ 

d^Cn-DyCn-l)  = g(n-l)  - I d^ (n-l)y (n-l+s) 

s=l 

Since  this  system  is  upper  triangular,  it  is  solvable  by  determining 
y(n-l)  from  the  last  equation,  y(n-2)  from  the  next-to-last  equation, 
and  so  on.  This  process  is  equivalent  to  developing  y (n-1) , y (n-2) , . . . , 
y(i)  from  (1.2.1)  by  backward  recurrence  starting  with  the  values  (1.2.3). 

The  noncriangular  cases  of  (1.2.4)  in  which  0 < j < ^ will  be  of 
greater  interest  to  us.  Our  primary  goal  is  to  find  a way  to  solve  these 
problems  by  performing,  in  effect,  a forward  recurrence  of  order  j 
followed  by  a backward  recurrence  of  order  £ - j . 

Definition  1.2.1  A nonsingular  boundary  value  problem  (1.2.4),  or  equi- 
valently (1.2.1)  to  (1.2.3),  having  initial  point  i , terminal  point  n, 
number  of  initial  conditions  j and  difference  operator 
D = [d^,dj^ , . . . ,d^]  , will  be  called  factorizable  provided  there 
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exist  difference  operators  A = [a^, . ,a^l  of  order  i and 

0 1 j 

B ^ = [b^  ^ » • • • ] of  order  k = I - j such  that 


(1.2.5) 


i,n-j-l  ^ i,n-j-l  i+j,n-l 


If  in  addition  the  sequences  and  free  of  zeros  and  the 

infinite  matrix  factorization 


(1.2.6) 


J J 


is  valid,  then  we  shall  say  that  the  difference  operator  D is  (i, j ) - 
factorizable . 

Several  remarks  may  be  made  here.  Note  that  ^ ^ and  are 

: J 

band  matrices  each  having  total  bandwidth  t + 1 , lower  bandwidth  j , 

and  upper  bandwidth  k . A nonsingular  and  factorizable  boundary  value 

problem  is  one  for  which  a finite  matrix  factorization  of  ^ ^ 

exists  such  that  the  left  factor  is  lower  triangular  with  total  bandwidth 

j + 1 and  the  right  factor  is  upper  triangular  with  total  bandwidth 

k + 1 . The  difference  operators  A^  and  may  be  arbitrary  as  long 

as  they  satisfy  the  condition  (1.2.5),  since  (1.2.5)  involves  only  a finite 

submatrix  of  each  of  A^  and  . Equation  (1.2.5)  is  illustrated 

in  Figure  4.  Since  ^ ^ is  nonsingular,  both  A^’'^  ^ and 

J J 

Bq"*"^  ’ ^ nonsingular.  Thus  the  entries  on  the  principal  diagonals 

of  A^’'^  ^ ^ and  ^ are  all  nonzero. 

J 0 

Noxa7  suppose  the  difference  operator  D € is  (i,j)  - factorizable; 

it  need  not  be  a nonsingular  operator.  Because  of  the  triangular  forms 
of  Aj  and  B ^ , equation  (1.2.6)  implies  equation  (1.2.5)  for  every 


n>i+j+l.  In  words,  every  leading  principal  submatrix  of  is 

j 

equal  to  the  product  of  the  corresponding  leading  principal  submatrices  of 
A^  and  . The  assumption  that  the  subsequences  a^  and 

jo 

are  free  of  zeros  is  equivalent  to  assuming  that  the  entries  on  the 


Figure  4.  The  factorization  (1.2.5). 
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principal  diagonals  of  the  infinite  matrices  and  are  all 


non- 


zero. It  follows  that  every  leading  principal  submatrix  of  is 


J 

nonsingular.  Thus  we  have  proved  the  following  theorem: 

Theorem  1.2.1  Let  D ^ be(i,j)-  factorizable . Then,  for  every 

n > i + j + max(j,£-j)  , the  boundary  value  problem  (1.2.4)  having  differ- 
ence operator  D,  initial  point  i,  number  of  initial  conditions  j and 

terminal  point  n is  nonsingular  and  factorizable.  j^j 

The  next  theorem  shows  it  is  possible  to  obtain  a solution  to  a fac- 
torizable boundary  value  problem,  in  theory  at  least,  by  performing  an 
appropriate  forward  recurrence  followed  by  an  appropriate  backward 
recurrence . 

Theorem  1.2.2  Let  (1.2.4)  be  a nonsingular  and  factorizable  boundary 

value  problem  with  initial  point  i,  number  of  initial  conditions  j,  term- 

i i 

inal  point  n and  difference  operator  D € . If  A and  B are 

difference  operators  of  order  j and  k = -£-j  , respectively,  such  that 
(1.2.5)  is  valid,  then  the  solution  of  (1.2.4)  is  identical  to  the  solu- 


tion of 


(1.2.7) 


i+j,n-l  i+j,n-l  _ i+j,n-l 

y ~ ^ 


0 


(A 


n-£,n-j-l  -1  n--£,n-j-l  n,n+k-l 


■)  D 


I 


y 


, i+j,n-l  , , ...  - 

where  z is  the  solution  of 


(1.2.8) 


i,n-j-l  i+j,n-l  _ i,n-j-l 
Aj  2 - g 


i.i+j-i  1,1+j-i 
0 


0 


are  nonsingular  because 


Proof  The  matrices  ^ ^ and  ^ 

— — ^ j 0 

i n*“  i ■“  2- 

D^’  is  nonsingular;  see  (1.2.5).  Therefore,  the  diagonal  elements 

.i»ti-i-l  , „i+i,n-l  , , . , . ->  . n 

or  and  (which  are  triangular  matrices)  are  nonzero 

and  (1.2.7)  and  (1.2.8)  are  nonsingular  linear  systems.  Furthermore, 
n-c.,  n-j-1 


A 


J 


is  nonsingular  because  it  is  a principal  submatrix  of 


A^’^  ^ . Therefore  the  right  side  of  (1.2.7)  is  unambiguously  defined. 

Now,  premultiplying  (1.2.7)  by  A^’'^  ^ ^ and  using  (1.2.5)  and 
(1.2.8),  we  obtain 


j^i,  n-j-1  i+j,n-l 

j ^ 


_ 

0 

0 

j 

, n-£, n-j-1  -1  n-£, n-j-1  n,n+k-l 

Consider  the  last  product  on  the  right  side  of  this  equation.  Its  block 
structure  is 


i,n-£-l 

j ; 0 

0 

^ : ^n-£, n-j-1 

‘ j 

^^n-£ , n- j -1  ^ ’ n- j -l^n , n-k-1 

^ J — 

where  0 in  the  matrix  denotes  a zero  submatrix  with  n - i - £ rows 
and  k(=  £-j)  columns,  0 in  the  vector  denotes  a zero  subvector  of 
length  n - i - £ , and  * in  the  matrix  denotes  a submatrix  which  needs 
no  further  identification  here.  Therefore,  the  product  in  question  is 


n-£, n-j-1  n,n+k-l 


equal  to  the  vector 


Substitution  into  the  equation  and  comparison  with  (1.2.4)  completes 
the  proof, 

An  important  remark  is  that  if  (1.2.6)  is  valid,  i.e.  if  D is  (i,j)- 
factorizable  with  left  factor  and  right  factor  , then 


n-^,n-j-l  _ n--C,n-j-l  n-k,n-l 

T .-u  • r,n-k,n-l  1 r / A ~1 

In  this  case  may  be  written  in  place  of  (A  ) 

j^n  t,n  i 1 (1.2.7).  The  truth  of  this  remark  is  apparent  from  Figure  4 

The  next  theorem  gives  necessary  and  sufficient  conditions  for  a 

difference  operator  to  be  (i,j)  - f actorizable . The  proof  also  shows 
that  (i,j)  - factorizations  may  be  constructed  by  Gaussian  elimination  in 
natural  order,  i.e.,  without  any  row  or  column  interchanges. 

Theorem  1.2.3  Let  D € . Then 

i)  A nonsingular  boundary  value  problem  having  initial  point  i, 

number  of  initial  conditions  j and  terminal  point  n is  factorizable  if 
every  leading  principal  submatrix  of  ^ ^ is  nonsingular. 

ii)  D is  (i,j)  - factorizable  if  every  leading  principal  submatrix 

of  is  nonsingular.  Furthermore,  this  condition  is  necessary  if  D 

is  nonsingular. 

Proof  Our  approach  is  to  apply  Gaussian  elimination  [10,  §2.1]  without  any 
row  or  column  interchanges.  Because  Gaussian  elimination  is  so  familiar, 
we  merely  sketch  the  proof.  The  proof  of  (i)  is  contained  in  the  first 
n - i - j stages  of  the  proof  of  (ii);  therefore  we  supply  only  the 
proof  of  (ii).  We  begin  by  proving  sufficiency  of  the  condition. 

The  first  row  of  is  used  to  annihilate  every  nonzero  element 

a 

below  the  first  element  in  the  first  colum.n.  This  is  possible  because 


the  pivot  element,  which  is  the  leading  principal  minor  of  order  one,  is 
nonzero.  Because  of  the  band  structure,  only  the  j rows  immediately  follow- 


ing  the  first  row  are  affected.  Also  (again  because  of  the  band  structure) 
the  trailing  coefficient  in  each  of  these  rows  is  not  changed.  The 
multipliers  used  in  these  annihilations  go  into  positions  2 through  j + 1 
of  the  first  column  of  the  lower  triangular  factor  . The  first  ele- 

ment is  unity  and  every  other  element  of  the  first  column  of  A^  is  zero. 

j 

This  completes  the  first  stage. 

Now  suppose  the  (r-l)st  stage  has  been  completed.  Let  the  notation 

6^  denote  the  transformation  of  the  original  resulting  from  the  first 

r - 1 stages.  The  leading  principal  submatrix  of  of  order  r is 

upper  triangular,  having  resulted  in  effect  from  the  corresponding  principal 

submatrix  of  the  original  by  Gaussian  elimination  in  natural  order. 

Its  diagonal  elements  are  all  nonzero,  by  hypothesis.  Therefore  the 

pivotal  element  is  nonzero  and  the  r-th  row  of  may  be  used  to 

annihilate  every  nonzero  element  below  the  r-th  element  in  the  r-th  column. 

Only  the  j rows  immediately  folloxd.ng  the  r-th  row  are  affected,  and  the 

trailing  coefficient  in  each  of  these  rows  is  not  changed.  The  multipliers 

used  in  this  annihilation  go  into  positions  r + 1 through  t + j of  the 

r-th  column  of  A^  ; unity  goes  into  the  r-th  position,  and  zeros  go  into 

every  other  position.  This  completes  the  r-th  stage. 

Let  be  the  infinite  upper  triangular  matrix  which  results  from 

the  full  sequence  of  transformations  on  . By  our  construction,  both 

and  B have  the  bandwidths  required  by  (i,j)  - f actorizability . 

Also  by  our  construction,  formally  multiplying  A^  on  the  right  by 

results  in  ; compare  (1.2.6). 

J 

In  order  to  prove  necessity  of  the  condition  in  part  (ii),  note  that 

the  factorization  produced  by  the  foregoing  process  is  such  chat  leading 

i i+i 

and  trailing  coefficients  of  both  A.  and  B 


are  nonzero  if  D is 
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nonsingular.  Therefore,  no  row  or  column  interchanges  could  be  made  with- 
out xd-dening  the  bandwidth.  Therefore,  Gaussian  elimination  in  natural 
order  is  the  only  possibility  for  arriving  at  the  factorization.  But 
Gaussian  elimination  in  natural  order  is  possible  only  if  every  leading 
principal  minor  is  nonzero.  This  completes  the  proof.  |_| 

Examination  of  the  proof  of  Theorem  1.2.3  shows  that  a unique  (i,j)- 
f actorization  does  not  exist,  since  a different  scaling  of  the  rows  is 
possible.  This  would  result  in  the  appearance  of  entries  on  the  principal 
diagonal  of  that  are  other  than  unity.  In  fact,  the  scaling  used 

in  the  (0, 1) -factorization  in  Giver's  second-order  algorithm  put  the 
sequence  (p ( 2) , p (3) , . . . ) on  the  principal  diagonal  of  the  right  factor; 
compare  §1.1.  Thus  the  principal  diagonal  of  the  left  factor  will  contain 


entries  other  than  unity,  in  general. 


1 . 3 A Preliiainarv  Development  of  the  Algorithm 


Let  D € . It  is  clear  that  an  extension  of  Diver's  algorithm 

for  computing  a particular  solution  y of  a given  difference  equation 
Dy  = g should  involve  the  posing  of  a boundary  value  problem  whose  solu- 
tion y^  is  an  approximation  of  y over  some  finite  subsequence  y^’^  . 
The  boundary  value  problem  will  have  some  number  j ^ 0 of  initial 
conditions  (depending  on  the  particular  solution  being  sought) . These 
will  consist  of  a subsequence  of  j consecutive  values  of  y . The 

terminal  conditions  of  the  boundary  value  problem  will  be  v (n)  = y (n+1) 

' n n 

= • • • = y (n+k-1)  = 0 , where  k = Z - j . 
n 

The  first  task  is  to  decide  how  to  choose  j . A necessary  pre- 
requisite is  a classification  of  possible  solution  types  when  the  operator 
is  of  order  higher  than  two.  Oliver  [19]  provided  an  early  extension  of 
Diver's  algorithm.  He  also  provided  the  following  definition  of  dominance. 
Suppose  we  are  given  two  sequences  a , b € S . Then  b dominates  a 
if  there  exists  some  i > 0 such  that 

(1.3.1)  |b(r+l)/b(r) I > | a (r+1) /a (r) | 

for  all  r > i . Oliver  considered  only  those  difference  operators  for 
which  a system  of  complementary  solutions  could  be  found  such  that  the 
complementary  solutions  are  linearly  ordered  under  this  definition  of 
dominance.  Furthermore,  he  considered  only  those  particular  solutions  of 
the  given  linear  difference  equation  which  could  be  compared  to  each  of 
the  complementary  solutions.  But  not  all  possible  particular  solutions 
can  be  so  compared;  for  example,  no  particular  solution  which  has  an 
infinite  subsequence  of  zeros  can  play  the  role  of  either  a or  b in 


(1.3.1). 
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In  §1.1,  we  adopted  a different  definition  of  dominance:  b 

dominates  a if 


(1.3.2) 


lim  a (r) 
b ( r ) 


This  differs  from  Oliver’s  definition  (which  is  equivalent  to  requiring 
that  ultimately  [a(r)/b(r)|  be  monotonically  decreasing).  Another 
difference  is  that  the  limit  of  |a(r)/b(r)|  under  Oliver's  definition 
may  be  greater  than  zero,  as  the  example  a(r)  = 1 + 1/r  and  b(r)  =1 
immediately  shows. 

Let  D € . The  set  /((D)  = {x€S|Dx=0}  will  be  called  the  kernel 

of  D . Obviously,  /((D)  is  a linear  subspace  of  5 . Let  us  suppose 
that  /((D)  has  dimension  t , so  that  there  exist  t linearly  indepen- 
dent sequences  . . . x^  ( /((D)  which  span  /((D)  . A sufficient 

condition  for  this  to  be  the  case  is  that  D be  nonsingular;  compare 
[L4,§12.1].  We  shall  say  that  the  set  x^,X2,...,x^  forms  a basis 
for  D . 

Definition  1.3.1  We  shall  say  that  D € is  totally  separable  (as 

r ^ «)  provided  there  exists  a basis  y^>y2>-*-»y^  D such  that 

y dominates  y in  the  sense  of  (1.3.2)  for  each  s , 1 < s < . 

s+1  s 

Any  basis  which  satisfies  the  conditions  of  Definition  1.3.1  will  be 
called  totally  ranked . Such  a basis  is  not  unique  because  y^  , for  example, 
could  be  added  to  each  of  y^j.-.y^  and  the  resulting  basis  would  also 
be  totally  ranked. 

For  each  s , 1 5 s < , let  K denote  the  linear  subspace  of  S 

s 

which  is  spanned  by  y]_’y2’'’’’^s  * particular,  = K(D)  . It 

seems  obvious,  and  it  will  be  proved  in  a more  general  setting  in  Chapter 
2,  that  every  totally  ranked  basis  generates  the  same  subspaces.  The 
subspaces  satisfy 


c ^2  c ...  c = K(D)  . 


Each  sec  inclusion  is  proper,  since  the  dimension  of  K is  s . This 

s 

ascending  chain  of  subspaces  induces  an  obvious  classification  of  the 
solutions  of  Dy  = 0 into  ^ distinct  types.  However,  we  are  interested 
in  classifying  the  solutions  of  Dy  = g for  nonzero,  and  indeed,  arbitrary 
g 6 S . Therefore,  we  frame  the  definition  so  as  to  cover  all  cases: 
Definition  1.3.2  Let  D € be  totally  separable  with  totally  ranked 
basis  • Then,  for  each  y € S , we  shall  say  that  the 

type  of  y (with  respect  to  D)  is  s , where  0 < s < t , and  we  shall 
write  type(y)  = s , provided  that 


(1.3.3) 


lim 


y(r) 


0 


and 


(1.3.4) 


lim  sup 


y(r) 

y^(r) 


> 0 


Furthermore,  if  (1.3.3)  is  true  for  s = 0 , we  define  type(y)  = 0 , 
and  if  (1.3.4)  is  true  for  s = I , we  define  type(y)  = Z . 

The  type  of  each  y ^ S is  the  same  regardless  of  which  totally 
ranked  basis  is  used.  This  will  be  proved  for  a more  general  classifica- 
tion of  sequences  in  Chapter  2.  Thus  we  are  justified  in  defining  the 
type  with  respect  to  D rather  than  with  respect  to  a particular  totally 
ranked  basis. 

The  conditions  for  Olver's  algorithm  [21]  suggest  that  Definitions 
1.3.1  and  1.3.2  are  appropriate  for  an  extension  of  the  algorithm  to 
arbitrary  order;  see  also  Cash  [ 6 ] . But  it  would  be  desirable  to  have  a 
classification  scheme  which  is  applicable,  for  example,  to  every  difference 
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equation  with  constant  coefficients.  Not  every  constant-coefficient 

operator  is  totally  separable,  however.  In  the  next  chapter  we  shall 

introduce  a more  general  class  of  difference  operators  (separable 

operators)  which  does  include  every  constant-coefficient  operator. 

Incidentally,  the  class  of  difference  equations  considered  by  Oliver  [19] 

includes  every  constant-coefficient  equation,  since  the  standard  basis 

(see  §2.1  below)  is  linearly  ordered  under  his  definition  of  dominance. 

Before  stating  the  next  theorem,  we  explain  what  we  mean  by  matrices 

and  determinants  associated  with  the  name  Casorati.  If  x.  ,x_,...,x  € S 

12  s 

are  arbitrary  sequences,  then  for  each  r > 0 the  s x s matrix 

X,  (r)  X (r)  

12  s 

X (r+1)  X (r+1)  X (r+1) 

12  s 

X (r+s-1)  X (r+s-1)  ....  x (r+s-1) 

Z s 

will  be  called  a Casorati  matrix.  Let  X = {x, ,x_,...,x  } . We  shall  use 
12  s 

the  concise  notation  [X] (r)  to  denote  the  Casorati  matrix  above.  Further- 
more, |x|(r)  will  denote  the  Casoratian  (of  X at  r) , defined  by 

|xi(r)  = det  [X] (r)  , r = 0,1,...  . 

If  D € is  nonsingular,  Casorati 's  theorem  [14, §12. 11]  states  that 

x^,x^,  ...,x^  forms  a basis  for  D if,  and  only  if,  |x|(r)  ^ 0 for 
all  r . 

Theorem  1.3.1  Let  y be  a solution  of  the  difference  equation  Dy  = g , 
where  D $ has  a totally  ranked  basis  Y = {y^ jy2> • • • • Let 


j = type(y)  . Then  y is  uniquely  determined  by  its  values  y(i),y(i+l), 
...,y(i+j-l)  provided  that  i > 0 is  a point  at  which  the  leading  prin- 

cipal minor  of  order  j of  |Yj(i)  is  nonzero. 

Proof  Assume  j = 0 ; then  the  condition  on  the  Casoratian  is  vacuous 
and  we  must  show  y is  uniquely  determined  without  knowledge  of  any  of 
its  values.  Suppose  z ^ y is  another  solution  of  type  zero;  then  there 
exist  unique  scalars  (not  all  zero)  such  that 

y(r)  = z(r)  + a^y^(r)  + •••  + a^y^(r)  , r = 0,1,,..  . 


Let  be  the  nonzero  scalar  of  highest  subscript.  Dividing  through  by 

y^(r)  and  allowing  r , we  see  that  = 0 , which  is  a contradiction. 

Next,  assume  j > 0 . Suppose  z ^ y is  another  solution  of  type  j 
such  chat 

z(i-l-r)  = y(i+r)  , r = 0,l,...,j-l  . 

There  exist  unique  scalars  . . . ,a^  6 F (not  all  zero)  such  that 

y(r)  = z(r)  + a^y^(r)  + ...  + a^y^(r)  , r = 0,1 


It  is  proved  easily  that  each  of 
to  that  of  the  previous  paragraph 
by  the  linear  system 


The 


..,a^  is 
remaining 


zero  by  a method  similar 
scalars  are  determined 


aiyi 


^ i » i+j -1 

+ a. y . 

J J 


0 


(where  0 stands  for  a column  vector  of  j zeros  and  we  have  used  the 
notation  developed  in  §1.2),  This  system  is  nonsingular  by  the  assumption 
on  the  Casoratian;  therefore,  = 0 , t = , which  is  a 

contradiction, 
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Theorem  1.3.1  confirms  the  adequacy  of  the  selection  of  j initial 

conditions  in  the  approximating  boundary  value  problem,  where  j = type(y) . 

In  ^3.3  we  will  demonstrate  that  stability  considerations  demand  exactly 

j initial  conditions,  in  general,  for  a solution  of  type  j. 

Now  let  us  suppose  that  a value  of  j has  been  selected,  where 

3 > type(y)  . Since  v = y and  y = 0 , the  linear 

' n n 

system  to  be  solved  is 


(1.3.5) 


^i,n-j-l  i+j ,n-l  ^ i,n-j-l 

j ^ 


i.l+j-l 

0 ^ 

0 


compare  (1.2.4).  The  value  of  n is  to  be  found  in  the  course  of  the 
computation.  In  view  of  Theorem  1.3.1,  y is  the  unique  solution  of  the 
infinite  linear  system 


(1.3.6) 


i+j  i 

= g - 


0 


We  assume  that  D is  (i, j )-factorizable,  i.e.,  that  linear  difference 

operators  $ P , and  . exist  such  that  . 

J 't-3  3 1 

As  shown  in  the  proof  of  Theorem  1.2.3,  successive  rows  of  this  factoriza- 
tion, and  indeed  of  the  triangularization  of  the  infinite  system  (1.3.6), 
may  be  generated  indefinitely  by  the  forw^ard  stage  of  Gaussian  elimination 
without  pivoting.  The  stability  of  this  process,  also,  will  be  discussed 
in  §3.3. 

In  consequence  of  Theorems  1.2.1  and  1.2.2  and  the  terminal  conditions 
^n,n+k  1 _ Q ^ solution  of  (1.3.5)  may  be  obtained  as  the  solution  of 


i+j ,n-l  i+j ,n-l 


^l+j,n-l 


(1.3.7) 


where  ^ is  a finite  subsequence  of  the  solution  of 


(1.3.8) 


. i i+j  i 


= a — 


i,i+j-i  i.i+j-i 

0 ' 

0 


Comparison  of  (1.3.8)  and  (1.3.6)  and  reference  to  the  identity 

= A show  that 
j J 


(1.3.9) 


T,  i"*"!  i+j 

B y = 2 


Thus  (1.3.7)  represents  the  finite  linear  system  that  results  from  (1.3.9) 
by  truncation  at  the  (n-i-j)th  row  and  column.  Since  and  z^"*"^ 

are  generated  by  Gaussian  elimination  without  pivoting,  the  solution  of 
(1.3.5)  for  any  n is  easily  programmed  for  an  automatic  digital  computer 
The  algorithm  comprises  two  stages:  ( i)  forward  elimination  (solution  of 

(1.3.8));  (ii)  back  substitution  (solution  of  (1.3.7)).  The  forward 
elimination  can  be  continued  arbitrarily  far.  The  back  substitution  is 
numerically  equivalent  to  back^^^ard  recurrence,  and  proceeds  ontil  we 
return  to  the  initial  point  i . 

The  question  of  the  stability  of  the  back  substitution  is  easily 

settled.  For  totally  separable  operators  we  have  the  following  theorem: 

Theorem  1.3.2  Let  D z have  a totally  ranked  basis  Y = {y^ , . . . . . y ^ 

Suppose  i > 0 and  j £ {0, 1, . . . ,£-l}  are  such  that  D is  (i,j)  - 

factorizable  and  the  leading  principal  minor  of  order  j of  |Y|(i)  is 

nonzero.  If  = A^’B^'*’^  is  an  (i,j)  - factorization  of  D , then 

J J 

is  totally  separable.  Furthermore,  if  z € S has  type  j or  less  with 

respect  to  D , then  z^"^^  has  type  zero  with  respect  to  . 

Proof  Suppose  j = 0 . Then  A^  is  a diagonal  matrix  with  nonzero  ele- 
ments on  the  diagonal;  see  Definition  1.2.1.  For  each  y ^Y  we  have 
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= 0 , i.e. 

I 

I a^(r)b  (r)y(r+s)  = 0 , r = i,i+l,...  . 
s=0 

Therefore, 

I 

I b (r)y(r-i-s)  = 0 , r = i,i+l,  . . . 
s=0  ® 

which  proves  {y^|y  ^ Y}  is  a totally  ranked  basis  for  . It  follows 

that  has  type  zero  with  respect  to  B^  whenever  z has  type  zero 

with  respect  to  D . 

Suppose  0<  j<  Z . Let  U = {y^>y2» • • • »Yj } • Since  |u|(i)  ^ 0 
it  is  possible  to  find  a linear  combination  x of  y^yV2’‘’'yy^  such 
that  x(i) ,x(i+l) , . . . ,x(i+j-l)  take  on  any  prescribed  set  of  values. 
Therefore,  appropriate  linear  combinations  x^  , s = j+1 , j + 2,  . . . , 

can  be  found  such  that  y (r)  + x (r)  vanish  for  r = i, i+1 , . , . . i+j -1  . 

^s-j  ” ^s  ’ ^ ~ j+1 , j + 2,  . . . . It  is  easy  to  verify  that 

{Yj^j  • • • jYj  > , . . . , v^}  is  a totally  ranked  basis  for  D . Now  let 
^ ^ > "'^2 ’ ' ’ ' ’ • In  view  of  Theorems  1.2.1  and  1.2.2,  the  remark 

following  Theorem  1.2.2,  and  the  facts  that  Dv  = 0 and  v(i+r)  = 0 
for  r = 0,1,...,  j-1  , we  conclude  that  ^ satisfies 


i+j,n-l  i+j,n-l 

Bq  V 


0 

„n-k,n-l  n,n+k-l 

\ 


for  every  n > i + j + max(j,k)  ; compare  equation  (1.2.7).  On  letting 
n ->  oo  this  proves  B^'*'^v^'''^  = 0 . Thus  {vj"'*”'^ , v^"^^  , . . . , v^"^^  } is  a 
totally  ranked  basis  for  B^~^^  . Finally,  let  z be  of  type  j or  less 
with  respect  to  D . Let  s £ {l,2,...,k}  . Then 


— V oo 


1 


Z (r)  ^ z(r) 

V (r)  y (r){l+o(l)} 
s j+s 


as 


i'^i  T>  ^”^3 

which  shows  chat  z has  type  zero  with  respect  to  3 . Lj  j 

I 

In  consequence  of  Theorem  1.3.2  and  our  remarks  in  §1.1,  when 
j > type(y)  the  back  substitution  stage  of  the  algorithm  is  stable,  at 
least  for  all  r in  a sufficiently  large  range. 

It  remains  to  determine  the  appropriate  value  of  the  truncation 
parameter  n . It  is  necessary  to  introduce  Green ' s formula  for  linear 
difference  operators;  see,  e.g..  Miller  [13].  We  apply  Green's  formula 
not  CO  the  operator  D but  to  Che  operator  . Our  development  is 

essentially  the  same  as  that  of  Cash  [ 6 ] . 

Let  B = [b-,b  ,...,b  ] i V . The  infinite  upper  triangular  matrix  | 

0 1 k k jj 

B = ^ which  is  obtained  from  B by  deleting  the  first  | 

k columns  and  transposing  the  result,  i.e.  whose  entries  are  defined  by 


(1.3.10)  b (r)  = b (r+s)  , 

s k-s 

/s 

is  a linear  difference  operator  of  order  k . We  shall  call  B the 
operator  which  is  adjoint  to  ^ , or  more  briefly,  the  adjoint  operator 
of  ^ . An  operator  and  its  adjoint  satisfy  a number  of  elementary 
properties  which  are  easily  verified  from  the  definition.  Among  these  are 


(1.3.11) 

bP  . 

(1.3.12) 

(1.3.13) 

k 

i 
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(1.3.14) 


_ .pP+k.q+k  Tr 
®k  - ’ 


valid  for  q 2:  p and  p = 0,1,2,...  . 

Now  take  any  u,v  € S and  B $ P , and  assume  k > 1 . Green's 

K. 

formula  is  the  identity 

q-1 

(1.3.15)  I {v(r)Bu(r)  - u(r+k) Bv(r) } 
r=p 

p,p+k-l  Tr  p,p+k-l  p,p+k-l 
= (^V  ) ^ ^ 

q,q+k-l  Tr  q,q+k-l  q,q+k-l 
— (^v  ) n u 

valid  for  q>p+k,  p>0  . Green's  fonnula  may  be  verified  by  express- 
ing the  left  side  of  (1.3.15)  in  the  form 


and  reducing  this  expression  with  the  aid  of  the  definition  (1.3.10)  and 
the  identities  (1.3.11)  - (1.3.14). 

Theorem  1.3.3  Let  be  an  (i , j ) -factorization  of  D ^ , 

J 1 "4 

where  D is  totally  separable  and  j < ^ . Let  y be  a solution  of 

Dy  = g , where  type(y)  5 1 and  g € S . Let  ^ and  be 

n 

the  corresponding  solutions  of  (1.3.5)  and  (1.3.8).  Also,  let  s be  a 

s4“X 

fixed  integer  satisfying  i+j  +k-l<s<n  and  w = (w  (s+1) , 

s s 

w (s+2),...)  be  the  solution  of 
s 


(1.3.16) 


;>s-k+l  s+1 
B,  w 
k s 


^s-k+i,s  s-k+1,3 
w 


where  k = -£  - j and 


(1.3.17) 


W,  = (0,0,. ..,0,1)  . 


Then 


(1.3.18) 


n-1 

y^(s)  = 5^  w^(r)z(r)  78^(3)  , 


where  ^^(s)  is  the  (s-i-j+l)st  leading  coefficient  of  . 

Remark:  Equations  (1.3.16),  (1.3.17)  are  equivalent  to  the  initial  value 

problem 


k 

I b (r)w  (r+t)  = 0 , r = s-k+1 , s-k+2 , . . . 
t=0  ® 


with  initial  conditions  w (s-k+1)  = • • • = w (s-1)  = 0 , w (s)  = 1 • 

s s s 

compare  (1. 2.1)-(1. 2.3)  and  (1.2.4a). 

Proof:  Let  p=s-k+l,q=n,v=w  and  u = y in  (1.3.15) . By 

s n 

virtue  of  (1.3.7),  (1.3.16)  and  (1.3.17),  the  left  side  of  (1.3.15) 
becomes 

n-1 

I w (r)z(r)  . 
s 

r=s 


. q,q+k-l  n,n+k-l  , , , . , - , ^ 

Since  u = V = 0 , the  second  term  on  the  right  side  of 

' n 


(1.3.15)  is  zero.  Thus,  expanding  the  right  side  and  using  (1.3.17), 


we  derive 
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, P+k-l^Tr^p , p+k-l^p , p+k-1 


(0,0,...,0,1) 

t>o(s-k+l)  ...  b^_^(s-k+l) 

y (s-k+1) 
n 

1 

O' 

o 
cn  ♦ 

1 

yjs) 

= bg(s)y^(s)  . 


Equating  the  left  and  right  sides  thus  arrived  at,  we  obtain  (1.3.18). 
In  §3.2  with  appropriate  conditions  it  will  be  proved  that 


lim  y (r)  = y(r) 
n-^  n 


for  each  fixed  r > i + j . Thus,  allowing  n « in  (1.3.18),  we  have 


(1.3.19) 


y(s)  = ^ w (r)z(r)/b  (s) 

s 0 

r=s 


for  each  s>i+j+k-l.  Similarly,  if 


Ti  (s ) = y(s ) - y (s ) 
n n 


is  the  truncation  error  incurred  in  accepting  the  approximation  y (s ) 

n 


for  y(s)  , we  have 


(1.3.20) 


T]  (s)  = I w (r)z(r)/b  (s) 


r=n 


for  each  s>i+j+k-l  . 

Let  us  compare  (1.3.20)  with  (1.1.2),  the  truncation  error  expansion 
for  the  original  algorithm  of  Olver.  Suppose  1=2  and  j = k = 1 . 


The  difference  operator  which  is  implicit  in  (1.1.1)  is  D = [dQ,d^,d2] 


where 


dg(r)  = a(r+l)  , d^(r)  = - b(r+l)  , d^(r)  = c(r+l) 

for  r = 0,1,2,...  . Similarly,  the  right  side  of  (1.1.1)  is,  in  our 
present  notation, 


g(r)  = d(r+l)  , r = 0,1.2,...  . 


The  rov;s  of  the  (0 , 1) -factorization  are  scaled  so  that 


bg(r)  = p(r+l)  , b^(r)  = - p(r) 


for  r = 1,2,...,  where  = [tig>b^]  and  p = (p  (0)  , p (1)  , . . . ) is  the 

sequence  defined  in  §1.1.  The  solution  z^  of  the  forward  elimination 
stage  of  the  algorithm  is  given  by 


z (r)  = e(r)  , r = 1 , 2,  . . . ; 

compare  (1.3.8).  Finally,  it  may  be  verified  readily  that  the  solution 

of  the  adjoint  equation  Bw  (r)  = 0 , r = s,s+l,...  such  that  w (s)  = 
-j  - s = 

is  given  by 


w (r) 
s 


p(s)p(s+l) 

p(r)p(r+l) 


s , s+1 , . . . 


Making  these  substitutions  into  (1.3.20)  and  exchanging  the  roles  of  r 
and  s , we  obtain  (1.1.2)  . 

The  simple  form  of  the  equation  for  w (r)  depends  solely  on  the 

s 

particular  scaling  of  the  ( 0, 1) -factor iza tion  used  by  Olver  and  the  fact 
that  k = 1 . A similar  simplification  can  be  obtained  for  difference 
equations  of  arbitrary  order  provided  only  that  k = 1 . It  stems  from 
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the  fact  that  the  order  of  the  adjoint  equation  is  one,  hence  all  solu- 
tions can  be  expressed  as  multiples  of  a single  particular  solution. 

These  conclusions  were  also  arrived  at  in  Cash’s  paper  [6]. 

The  expansions  (1.3.19)  and  (1.3.20)  may  be  employed  in  the  estimation 
of  the  optimal  value  of  n such  that  a given  termination  criterion  is  met. 
Typical  termination  criteria  would  be  to  require  one  or  more  values  of 
y^(r)  in  the  range  i+j<r5m  to  approximate  y(r)  to  a specified 
absolute  or  relative  precision.  Since  the  back  substitition  stage  is 
stable,  controlling  the  absolute  or  relative  error  only  at  the  point  r = m 
is  often  sufficient.  In  practice  equations  (1.3.19)  and  (1.3.20)  are 
employed  by  replacing  the  upper  limit  by  n + v , where  v ^ 1 is  a 
parameter  chosen  to  suit  the  problem  at  hand. 

The  solution  of  (1.3.16),  (1.3.17)  is  obtained  by  for\^7ard  recurrence. 
The  stability  of  this  process  will  depend  on  the  behavior  of  for 

large  r relative  to  the  other  solutions  of  the  adjoint  equation.  In 

general  there  is  no  reason  to  expect  that  w (r)  will  not  contain  a com- 

s 

ponent  of  a maximal  solution.  Thus  the  solution  of  (1.3.16)  will  normally 
be  stable,  although  there  could  be  an  initial  range  where  some  unstable 
rounding  error  propagation  takes  place. 


CHAPTER  2.  CLASSIFICATION  OF  SEQUENCES 


2.1  Constant-Coefficient  Operators 

Let  D = . . .d^]  - be  a constant-coefficient  operator.  By 

this  we  mean  that  there  exist  6 ,6  ,...,S„  ^ F such  that  6 0 , 

5,  7^  0 and  d (r)  = 6 for  all  r . We  assume  without  loss  of 
■L  s s 

generality  that  5^  = 1 in  every  constant-coefficient  operator.  Since 
every  constant-coefficient  operator  is  nonsingular,  by  our  definitions,  the 
kernel  K(D)  has  dimension  Z ; see  discussion  preceding  Definition  1.3.1. 

A particular  basis,  which  we  shall  call  the  standard  basis , exists 
for  each  constant-coefficient  operator;  see  [ 11  ,§§1 3 .0-13 . 1]  . The  polynomial 

I 1-1 

P(C)  = r + + •••  + 

is  known  as  the  characteristic  polynomial  for  D . The  standard  basis  is 
the  set  of  sequences 

1 = {(r°^\^)  € S I P(\)  =0  and  m 6 {0, 1,  . . . ,M.(X)-1}} 

where  M-(X)  denotes  the  multiplicity  of  the  root  X . We  are  using  the 
notation  (r'^^X^)  here  as  an  alternative  way  of  denoting  the  sequence 

(p  ,X, 2 X ,3  X , , . . ) 
m 

where  p„  = 1 and  p =0  when  m > 0 . 

0 ra 

The  standard  basis  of  a constant-coefficient  operator  is  not  necessarily 
totally  ranked;  recall  Definition  1.3.1.  For  example,  if  X^  and  X^ 
are  two  distinct  zeros  of  the  characteristic  polynomial  such  that 
|X^!  = jX^I  then  |X^|  = | X^ | for  all  r . Neither  solution  dominates 


the  other. 
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Let  us  x^rite  a < b whenever  a,b  c S are  such  that  a is  dominat- 
ed by  b . We  shall  say  that  a and  b are  separated  (at  infinity)  if 
a<b  or  b<a.  A totally  ranked  basis  is  linearly  ordered  by  < 

but  in  general  the  standard  basis  is  not,  as  we  have  just  seen.  However, 
we  can  introduce  a linearly  ordered  partition  of  the  standard  basis . We 
shall  proceed  to  do  this. 

m,  m„ 

1 r 2 r 

It  is  apparent  that  (r  < (r  if,  and  only  if,  either  (i) 

|l-ll  < (ii)  = 1^2!  ™1  ^ denote  the 

number  of  roots  of  P(^)  =0  of  distinct  absolute  value,  and  let 

X , X , ...,\  be  representative  roots  of  P(^)  = 0 such  that  |X  | < |k  j 

1 2 O’  s s+ i 

for  each  s . Define  u.  = max{!-i.(X)!  |Xl  = !X  1}  . We  define  the  sets 

' s s 

= {(r'^x’')  t a|  |Xl  = |Xj}  , s = l,2,...,c 

and,  for  each  s , the  sets 

A = {(r^X^)  € A |p.(X)  > m}  , m = 0,1,  . . . ,p  -1  . 

S 5 Til  S 3 

It  may  be  verified  readilv  that  the  sets  A are  disjoint  and  nonempty 

s , m 

and 

u -1 
O’  ' s 

(2.1.1)  A = U U A 

1 n s,m 

s=l  m=0 

The  sets  A are  linearly  ordered  according  to  the  definition: 

S ) HI 

XI  < A , , if,  and  only  if,  either  (i)  s < s’  or  (ii)  s = s'  and 

s , m s , m 

m < m ' . 

An  illustrative  example  is  afforded  by  the  constant-coefficient  equa- 


tion 


y(r+8)  - (4+4V'2)y(r+7)  + (15+16v'I)y (r+6)  - (48+12V2)y (r+5) 

- 12y(r+4)  + (192+64v'2)y (r+3)  - ( 208+256Vl)y (r+2)  + ( 256+19 2V2)y (r+1) 

- 192y(r)  = 0 . 

Its  characteristic  polynomial,  in  real  factored  form,  is 
(2.1.2)  (4-1)(C-2)(C+2)(C^-2\/2?+4)^(5-3)  = 0 . 

The  roots  of  (2.1.2)  include  the  complex  numbers  V2(l+i)  and  V2(l-i)  , 
each  of  multiplicity  two.  The  subsets  of  the  linearly  ordered  partition 
in  this  example  are 


‘1,0 


A 


2,0 


^2,1 


3,0 


= {2^(-2)^V2U-i)^V2"(l+i)^} 


{rV2'"(l-i)’',rV^’'(l+i)’'} 


{3^} 


Note  that  if  a 6 ^ and  ^^^"2  1 ia(r)/b(r)|  = r 


-1 


V7hen  r > 1 


In  contrast,  if  a 6 A and  b 6 A , where  s < s’  then  !a(r)/b(r)| 

s s 

contains  a decreasing  exponential  factor. 

In  general,  if  a,b  6 A are  such  that  a < b , then  the  separation 
ratios  ja(r)/b(r)|  decrease  as  the  product  of  a power  of  r ^ and  an 
exponential  function  in  r . If  the  exponential  function  is  identically 
equal  to  one,  as  will  be  the  case  whenever  a and  b are  standard  basis 
solutions  which  correspond  to  characteristic  roots  of  equal  magnitude, 
then  we  shall  say  that  a and  b are  algebraically  separated . Otherwise, 
we  shall  say  that  a and  b are  exponentially  separated . 

Figure  5 is  obtained  by  plotting  the  roots  of  (2.1.2)  in  the  complex 


plane  and  drawing  circles  through  them,  centered  at  the  origin.  There  are 
three  circles,  corresponding  to  the  distinct  absolute  values  of  the  roots. 
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The  double  roots  are  indicated  on  the  circle  of  radius  two;  it  is  help- 
ful to  imagine  them  duplicated  on  a second  circle  of  radius  two  that  is 
raised  slightly  out  of  the  plane.  Let  us  number  the  four  circles,  thus 
arrived  at,  hy  increasing  three-dimensional  distance  from  the  origin. 

This  determines  geometrically  the  proper  linear  ordering  of  the  four  sub- 
sets A of  the  partition.  The  ordering  of  elements  within  any  one  of 
s > m 

the  subsets  is  unimportant  for  our  purposes. 

Clearly,  this  geometric  indication  of  the  linear  ordering  of  the  sub- 
sets of  the  partition  can  be  extended  to  any  constant-coefficient  operator. 


Figure  5.  The  roots  of  (2.1.2).  The  double  roots  ''/2(l±i)  are 


denoted  X,  \ . 


.1  Let  (2.1.1)  denote  the  linearly  ordered  partition  of  the 


. r.e  ere-  _ . 


standard  basis  for  a given  constant-coefficient  operator.  Then: 


i)  If  a ^ and  b £ '•  , , , then  a < b if,  and  only  if, 

^ s , n 


ii)  If 

iii)  If 


a,b  c ' , then  a(r) ' = ib(r)-  for  all  r . 

3 ,ni 


a,  , a , . . . , a , b € A 
1 2 q ; 


where  the  a's  are  distinct,  and 


7 ’ ‘ 


..,c  € F are  all  nonzero,  then 


(2.1.3) 


lin  sun  T t a (r)/b(r)  > 0 . 

P P 

r-^  p=l 


Proof:  Parts  (i)  and  ( ii)  are  trivial  consequences  of  our  definitions. 

Turning  to  oart  (iii),  we  note  that  the  sequence  c a (r)/b(r)  is 

p=l  p p 

not  identically  zero  (since  the  a's  are  linearly  independent).  Let 

,\  be  the  c'naracteristic  roots  associated  with  the  set  A 
1 ^ q s,E 

such  that 


O D 


p = 1 , 2 , . . . , q 


n,  r 


and  b(r)  = r X , Since  the  X's  lie  on  the  same  circle  in  the  complex 
plane,  there  exist  real  constants  lying  in  the  half- 
open real  interval  [-"',“)  such  that 


a^(r)  = , p = 1,2, ...,q 


and  b(r)  = r^  \,"^exp(irp)  , i - -1  . Thus 


I t a (r)/b(r)  = i z 
p=l  P P p=l  P 


where  ' = ';  - 3 • The  v's  are  distinct  real  numbers  (modulo  2n) 

P P 
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because  the  sequences  a are  assumed  distinct. 

1 Z q 

The  right  member  of  the  last  equation  represents  a so-called  almost 
periodic  function  of  r when  r is  regarded  as  a continuous  real  var- 
iable. To  complete  the  proof  we  appeal  to  the  following  lemma,  the  proof 
of  which  is  given  in  the  Appendix. 

Lemma  2.1.1  Let  a,  be  nonzero  complex  numbers  and  r, 

1 2 q 12 

distinct  real  numbers  (modulo  2tt)  . Define  the  infinite  sequence 


ir"v* 

(2.1.4)  x(r)  = y a e P , r = 0,1,2,...  . 

1 P 

p=l 


Then  lim  sup^_^^j x(r)  | > 0 . 


An  equivalent  statement  to  inequality  (2.1.3)  is  that  there  exists  an 
infinite  subsequence  of  the  sequence 


q 

I a a (r)/b(r)  , r = 0,1,2, ..  . 

p=l  P P 


which  is  bounded  away  from  zero.  Intuitively,  this  means  that  it  is 

impossible  to  create  a new  sequence  with  qualitatively  different  asymptotic 

behavior  merely  by  forming  a linear  combination  of  standard  basis  solutions 

of  similar  asymptotic  behavior.  Note,  however,  that  unlike  the  original 

sequences  infinitely  many  zeros  may  occur  in  the  new  sequence.  For  example, 
r r 

if  2 and  (-2)  are  sequences  in  a standard  basis,  then  every  other 

r r 

term  of  the  sequence  2 + (-2)  is  zero. 


We  relegate  the  proof  of  Lemma  2.1.1  to  the  Appendix  because  the  theory 
of  almost  periodic  functions  is  not  germane  to  the  rest  of  this  thesis. 


2.2  Separable  Operators 


In  this  section  we  introduce  a subclass  of  the  general  set  of  infinite 


detail  the  general  structure  of  the  kernel  of  an  operator  in  this  sub- 
class. This  will  prepare  the  way  for  the  presentation,  in  the  next 
section,  of  a general  classification  of  all  sequences  in  S with  respect 
to  a given  linear  difference  operator.  In  Theorem  2.1.1  we  presented 
three  properties  of  the  standard  basis  of  a constant-coefficient  operator. 
These  properties  suggest  the  following  extension  of  Definition  1.3.1: 
Definition  2.2.1  Let  D i V ^ be  such  that  fC(D)  has  dimension  t . 

Then  D will  be  said  to  be  separable  if  there  exists  a basis  X for  D 
which  can  be  partitioned  into  nonempty  disjoint  subsets  in 

such  a way  that  the  following  three  conditions  are  all  satisfied: 

i)  If  X € X and  y € X , where  s < t , then  lim  x(r)/y(r)  = 0 

S t r^ 

ii)  If  x,y  6 X for  some  s , then  0 < lim  inf  lx(r)/y(r)| 

s r->" 

and  lim  sup^_^^|  x (r) /y  (r)  | < °°  ; 

iii)  If  X is  any  linear  combination  of  sequences  from  X^  , other 
than  X = 0 , and  y £ X^  , then  lim  sup^^_^ ! x(r) /y  (r)  | > 0 . 

Let  D be  a separable  operator.  A basis  X for  D which  is  such 

that 


where  the  X satisfy  the  conditions  of  Definition  2.2.1  will  be  said  to 
s 

be  optimally  ranked . Let 


CT 


(2.2.1) 


X = U X , 

1 ^ 
s=l 


(2.2.2) 


where  the  subscript  n is  positive  (because  X 

s s 


is  nonempty) . 
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Obviously,  we  have 
(2.2.3) 


t-  I 


s=l 


n 

s 


In  Lemma  2.2.1  below  we  will  show  that  two  distinct  optimally  ranked  bases 
for  a given  separable  operator  have  the  same  number  (cr)  of  disjoint  subsets 
and  the  same  number  (n^)  of  sequences  in  corresponding  subsets. 

An  optimally  ranked  basis  exhibits  the  full  range  of  possible 
asymptotic  behavior  of  solutions  to  the  corresponding  homogeneous  differ- 
ence equation.  Conditions  (i)  and  (ii)  of  Definition  2.2.1  establish  the 
ranking  of  sequences  within  the  basis.  If  (2.2.1)  is  a particular  optimally 
ranked  basis,  we  shall  say  that  x € X is  a sequence  of  type  s provided 
that  X € . More  generally,  we  shall  say  that  x € K(D)  is  a sequence 

of  type  s orovided  that  x 0 and  there  exist  scalars  a € F such  that 

^ ‘ p»q 

n 


(2.2.4) 


X 


s 


I 

P=1 


a X 

p.q  p>q 


where  a ^0  for  some  q , 1 5 q 5 n . Such  a representation  of  x 
s,q  ^ ^ s 

exists  and  is  unique  for  every  nonzero  x € K(D)  . Again,  in  consequence 
of  Lemma  2.2.1  below,  if  we  change  to  a different  optimally  ranked  basis, 
the  type  of  every  sequence  in  /((D)  remains  invariant.  Condition  (iii) 
of  Definition  2.2.1  presages  this;  it  requires  that  no  nontrivial  linear 
combination  of  basis  sequences  of  the  same  type  may  be  dominated  by  a 
basis  sequence  of  that  type. 

Let  L(Y)  denote  the  linear  subspace  of  S which  is  spanned  by  Y 

where  Y c S may  be  any  finite  subset  of  sequences.  If  (2.2.1)  is  an 

optimally  ranked  basis  for  a separable  operator  D,  then  we  shall  say  that 
s 

the  subspace  I ( U X ) is  the  subdominant  subspace  of  type  s for  D . 

- p 

p=i 


These  subspaces  are  independent  of  the  basis;  again,  see  Lemma  2.2.1 

s 

below.  If  i denotes  the  dimension  of  L(  U X ) , then 

" P=i  p ■ 


(2.2.5) 


' I 

p=l 


compare  (2.2.2)  and  (2.2.3).  In  addition,  we  shall  say  that  {0}  , the 
linear  subspace  of  dimension  zero  (whose  only  sequence  is  the  zero  sequence) 
is  the  subdominant  subspace  of  type  0. 

Let  U be  one  of  the  subdominant  subspaces  for  D . The  complementary 
subspace  of  (J  in  f((D)  , i.e.,  the  subspace  1/  of  K(D)  such  that 


(2.2.6) 


K(D)  = a © 1/  , 


will  be  called  the  corresponding  dominant  subspace.  The  symbol  ® on  the 
right  side  of  (2.2.6)  means  that  K(D)  is  the  direct  sum  of  the  subspaces 

U and  1/  ; that  is,  every  nonzero  sequence  x 6 /((D)  is  uniquely 

expressible  in  one  of  the  forms  x = u , or  x = v , or  x = u+v  , 

where  u € U and  v 6 1/  . Clearly,  if  U is  the  subdominant  subspace 

a 

of  type  s,  0 5 s < a , then  1/  = L(  U X ) ; and  if  U is  the  sub- 

p=s+l  P 

dominant  subspace  of  type  a,  then  1/  = {0}  . 

Before  stating  and  proving  Lemm.a  2.2.1,  we  note  that 

s 

L(  U X ) = L(X. ) © L(X  ) a • • • ® L(X  ) 
p=i  p ^ ^ 

Thus,  an  alternative  characterization  of  a sequence  of  type  s in  f((D)  is 
the  following:  x € K(D)  is  of  type  s,  1 5 s < a , provided  that 


s 


r 


X 

p 


6 L(X  ) 
P 


X 0 ; 

s 


compare  (2.2.4) . 
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Lemma  2.2.1 

Let 

and 

CT  O’ 

X = U X = U {x  ,x  ^ , . . . ,x  } 

, s , s,l  s,2  s,m 

S=1  s=l  s 

T T 

Y = U Y = U {y  ,y  , . . . ,y  } 

.s  ,s,ls,2  s,n 

S=1  s=l  s 

be  distinct  optimally  ranked  bases  for  a separable  operator  D ^ . 


Then  cr  = t 

and  m = n for  all  s . Furthermore, 

s s 

s s 

L(UX)=L(UY),  s = 1,2,  . . .,  <j  , 
, s , s 

P=1  P=1 

and  the  type  of  every  nonzero  sequence  z € K(D)  is  the  same  with  respect 
to  X as  it  is  with  respect  to  Y , 

Proof  Let  u S and  v € X^  , so  that  relative  to  X the  types  of  u 


and  V are 

s and  t . Let  the  types  of  u and  v relative  to  Y 

be  s'  and 

t'  . Then  there  exist  unique  representations 

and 

s ' 

u = } u.  , u.€L(Y.),  u,^0 

jii  j : j 

t' 

V = 1 V.  , V.  6 L(Y.)  , V , 7*  0 . 

jil  3 J : 

Furthermore,  by  condition  (i)  of  Definition  2.2.1,  we  have,  as  r , 

(2.2.7)  u(r)  = u^,(r)  + o{y(r)>  , y ^ , 


and 


(2.2.8) 


v(r)  = v^,(r)  + o{y(r)}  , y 6 , . 

i)  First  we  prove  that  s'=t'  if  s=t,i.e.,if  u and  v are  of 
the  same  type  with  respect  to  X , then  they  are  also  of  the  same  type 
with  respect  to  Y . 

Since  s = t , the  sequence  u(r)/v(r)  is  asymptotically  bounded 
away  from  zero;  compare  part  (ii)  of  Definition  2.2.1.  Suppose  that 
s'  < t'  . We  shall  obtain  a contradiction  by  showing  that,  under  this 
assumption,  u(r)/v(r)  has  an  infinite  subsequence  which  converges  to 
zero.  Choose  any  y ^ Y^ , . Then  we  have 

lim  u(r)  ^ Q 
r-x»  y(r) 


From  (2.2.8),  we  have 


v(r) 

y(r) 


v^,(r) 

y(r) 


+ 0(1) 


and,  using  the  triangle  inequality,  an  elementary  property  of  the  limit 
superior,  and  condition  (iii)  of  Definition  2.2.1,  we  obtain 

lim  sup  j V (r ) 

|y(r) 

Thus  there  exists  an  infinite  subsequence  of  v(r)/y(r)  that  is  bounded 
away  from  zero,  and  therefore  also  a subsequence  of 


lim  sup 

, Vi  y 

r^ 

y(r) 

u(r)  ^ u(r) /y(r) 
v(r)  v(r)/y(r) 


that  converges  to  zero.  This  supplies  the  needed  contradiction. 
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A similar  contradiction  is  obtained  if  we  suppose  that  s'  > t'  . Hence 
s'  = t'  . 

ii)  Next  we  prove  that  s'  < t'  if  s < t . 

Since  s < t , the  sequence  u(r)/v(r)  converges  to  zero.  Suppose 
that  s'  > t'  . We  shall  derive  a contradiction  by  showing  that,  under 
this  assumption,  u(r)/v(r)  has  a subsequence  which  does  not  converge  to 
zero.  Choose  any  y € , . First,  by  using  (2.2.7)  and  an  argument 

similar  to  that  used  in  part  (i)  of  this  proof,  we  have 


lim  sup 
>00 


u(r) 

y(r) 


> 0 . 


Thus  there  exists  an  infinite  subsequence  of  u(r)/y(r)  that  is  bounded 
away  from  zero. 

Assume  now  that  s'  > t'  . Then 

lim  v(r)  ^ Q 
r-»oo  y(r) 

In  this  case  it  follows  that  there  exists  a subsequence  of 


u(r)  _ u(r) /y (r) 
v(r)  v(r)/y(r) 


that  diverges  to  « , which  supplies  the  needed  contradiction. 

Alternatively,  assume  that  s'  = t'  . Using  (2.2.8),  we  have  for 
sufficiently  large  r 


v(r) 

Vj.,  (r) 

y(r) 

y(r) 

Furthermore,  using  condition  (ii)  of  Definition  2.2.1,  we  see  that  there 

exist  r >0  and  S > 0 such  that 
o 


,v^,(r) . 

'_J:: i <:  ? r > r 

' y(r)  j ' ^ ’ O 


Tnerefore,  for  sufficiently  large  r , we  have 


|V(r) 

'y(r) 


? + 1 


and 


iu(r) 

v(r) 


- (1+p) 


-l|u(r) 

'y(r) 


It  follows  that  u(r)/v(r)  has  a subsequence  that  is  bounded  away  from 
zero,  which  again  supplies  the  needed  contradiction. 

iii)  Next,  we  prove  that  i = z and  m = n for  all  s , and  also 

s s 

s s 

that  L(  j X ) = L(  U Y ) for  all  s . 
p-1  p=l 

Part  (i)  of  this  proof  implies  that  for  each  s , 1 5 s < a , the  type 
relative  to  Y of  every  sequence  in  is  the  same,  say  t^  . Using 

part  (ii),  we  have 


(2.2.9)  1 5 t < t < ••'  < t 5 T . 

12  O’ 

Furthermore,  that  s < t for  all  s , 1 5 s 5 o , is  easily  proved  bv 

s 

induction.  In  particular,  o < t and,  using  (2.2.9),  we  see  that  o 5 T 

o 

Inverting  the  roles  of  X and  Y and  repeating  the  argument,  we  see 
that  T S T . Thus  o = t , and  (2.2.9)  implies  t^  = s for  all  s . 

For  each  s , it  follows  that 


L(X,J---UXJ  c L(Y  IJ---UY  ) , 
IS  is 


and  by  symmetry. 


Thus 


L(X,U*  • -U'X  ) = L(Y.U-  • -UY  ) 

1 s 1 s 

for  each  s . In  particular,  l(X^)  = L(Y^)  . Since  X^  and  Y^  are 
bases  of  the  same  linear  space,  we  have  ~ easy  induction 

completes  the  proof  of  this  part. 

iv)  To  complete  the  proof,  let  z 6 K(D)  be  a nonzero  sequence  such  that 


z 

= 1 ^ 

u r 0 , 

u 

€ i.(X  ) , 

p=l 

P 

s 

P 

P 

and 

t 

z 

= 1 ^ 

V # 0 , 

V 

6 L(Y  ) . 

p=l 

P ’ 

t 

? 

P 

Assume  s < t , and 

let  V 

€ Y 

. Then 

t 

(-* 

< 

rt 

lim 

s 

V 

u ( r ) -V 
P ? 

(r) 

t-1  V (r) 

— n 

r-*-»  y(r) 

L 

_p=i 

y(r) 

p=Li 

— U . 

But  this  result  contradicts 

condition  (iii) 

of  Definition  2.2.1.  Hence 

must  have  s > t . 

A similar  contradict 

ion 

is  obtained  if 

we  assume 

s > t . Therefore 

s = t 

and 

the  proof 

of 

Lemma  2.2.1  is 

complete.  !_ 

In  general,  the 

actual 

determination 

of 

optimally  ranked  bases  for 

variable  coefficient 

difference 

operators  is 

a difficult  problem.  Some 

results  have  been  given  by  Wimp  [26].  This  reference  extends  earlier 
work  by  Birkhoff  and  Trjitzinsky  [2,3]  on  the  analytic  theory  of  singular 
difference  equations  whose  coefficients  possess  asymptotic  expansions 


of  a orescribed  form. 


All  equations  with  coefficients  rational 

in  r are  included,  for  example.  Wimp's  analysis  provides  a 

means  of  constructing  a basis,  which  he  calls  a canonical  set: 

see  [26,  Theorem  3.3  and  Definition  3.6].  The  canonical  set  is 

analogous  to  our  optimally  ranked  basis.  More  specifically,  there 

exists  a basis  {Yt  jY-j  • • • >7;?}  such  that,  for  each  s = , 

J-  z ^ 

y (r)  = c M (r){l+o(l)}  as  r ^ “ 

■ s s s 


where  c 7^  0 and 
s 


Q (r)  6 

M (r)  = e r ^(£n  r) 
s 


Here  9 is  a comnlex  number  and  p is  a positive  integer.  Q is  of 
s ‘ s s 

the  foirm 


Q^(r)  = a^r  In  r + u^r  + + •••  + 


where  0 is  an  integer,  p > 1 , and  p ,u..,  , . . . are  complex  numbers. 

U -L  p 

The  rich  variety  of  possible  separation  ratios  of  solutions  as  r -»■  «>  , 
compared  with  the  algebraic  or  exponential  separation  of  solutions  in  the 


constant-coefficient  case,  described  in  §2.1  above,  is  evident. 
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2.3  The  General  Classification 

In  §2.2  we  showed  that  every  nonzero  sequence  in  the  kernel  of  a 
separable  operator  has  a unique  t^/pe  with  respect  to  that  particular 
operator.  In  this  section  we  extend  this  classification  to  every  nonzero 
sequence  in  S . The  following  definition  should  be  compared  with 
Definition  1.3.2. 

Definition  2.3.1  Let  D € be  separable  with  optimally  ranked  basis 
a 

X = U X . Then,  for  each  y € 5 , we  sav  that  tvpe  (y)  = s , or  more 

3 = 1 ^ 

fully,  the  type  of  y with  respect  to  ^ is  s , provided  that  0 < s < c 
and  the  following  two  conditions  are  satisfied: 


i) 


If  X € X 


s+1 


then 


lim  y (r) 
r-*«>  x(r) 


0 ; 


ii) 


If  X T X then  = > 


lim  sup  I y (r) 

r-*oo  * X ( r ) 


> 0 


If  (i)  is  satisfied  for  s = 0 , then  (ii)  does  not  apply,  and  we  define 
type(y)  = 0 . Similarly,  if  (ii)  is  satisfied  for  s = cr  , (i)  does  not 
apply,  and  we  define  type(y)  = ct  . 

The  classification  of  sequences  introduced  in  Definition  2.3.1  allows 
us  to  treat  both  homogeneous  and  inhomogeneous  difference  equations  in  a 
uniform  manner  when  the  difference  operator  is  separable.  The  kernel  of 
the  operator  contains  every  solution  of  the  homogeneous  equation,  by  defini- 
tion. But  solutions  of  the  inhomogeneous  equation  with  arbitary  right-hand 
side  may  lie  anywhere  in  S . This  is  why  a classification  of  all  sequences 
in  S is  desirable. 

Theorem  2.3.1  Every  sequence  in  5 possesses  an  unambiguous  type  with 
respect  to  each  separable  operator  in  V ^ . 


Proof  (i)  First  we  show  that  the  type  with  respect  to  ^ fixed  optimally 


ranked  basis  is  determined  uniquely  by  Definition  2.3.1  for  every  nonzero 


sequence  in  S . 

Let  z € S . Assuming  that 


(2.3.1) 

lim  sup 

z(r) 

x(r) 

> 0 

for  at  least  one  sequence 

X ^ 

X , 

let 

X € X 

type  such  that  (2.3.1)  is 

satisfied.  Let  s 

is  satisfied  for  every  x 

6 X 

s 

, since 

z(r) 

z(r) 

x(r) 

x(r) 

x(r) 

x(r) 

be  a sequence  of  largest 
type  (ic)  . Then  (2.3.1) 


and  lim  inf ^ ^^jx(r) /x(r) | >0  as  a consequence  of  condition  (ii)  of 
Definition  2.2.1. 

If  s = CT  , then  type(z)  = a according  to  Definition  2.3.1.  There- 
fore suppose  that  s < ct  . By  the  maximality  of  x we  see  that 
lim  sup^_^  j z (r) /x (r) I = 0 whenever  type  (x)  > s or,  equivalently, 

(2.3.2)  = 0 whenever  type  (x)  > s . 

r-^  X ( r ) 

Accordingly,  type  (z)  = s . 

To  conclude  the  proof  of  part  (i),  let  us  suppose  that  (2.3.1)  is 

not  satisfied  for  any  x ^ X . Then  (2.3.2)  is  satisfied  for  every  x 6 X 

and  in  particular  for  every  x € X^  . Then  type  (z)  = 0 . 

a a 

(ii)  In  this  part  we  prove  that  if  X = U X and  X'  = U X[ 

s=l  ^ s=l  ^ 

are  two  distinct  optimally  ranked  bases  for  D then  the  type  of  every 

sequence  in  S is  the  same  with  respect  to  X and  X'. 

Choose  any  s , 1 < s 5 a , and  choose  x ^ X , x'  <:  X'  . There 

s s 


exist  X 

P 

Furthermore 


X , 1 < p 5 s , such  that 
P 

we  have 


x,+*--+x  ; see  Lemma  2.2.1 

i s 


x^(r) 

x(r) 


= 0(1)  + 


x (r) 
s 

x(r) 


as  r ^ . 


Since  x^  is  a linear  combination  of  sequences  of  type  s,  a brief  compute 
tion  and  application  of  condition  (ii)  of  Definition  2.2.1  yields 


lim  sup 


X (r) 

s 

x(r) 


< 


oo 


It  follows  that 


(2.3.3) 


lim  sup  X ' (r) 
r-»-“  x(r) 


and 


lim  inf 


xlrl 


0 whenever  type  (x')  = type  (x) 


x’(r) 

Since  the  roles  of  x and  x'  are  interchangeable,  we  also  have 

x(r) 


(2.3.4) 


lim  sup 


and 


< 


x'  (r) 

1 im  inf 


x*(rl 

x(r) 


> 0 whenever  type  (x')  = type  (x) 


Let  z € S and  suppose  type  (z)  = 0 with  respect  to  X . Thus 

z(r)/x(r)  -*•  0 for  every  x ^ X . Choose  x'  € X'  . Let  s = type  (x') 

and  choose  x C X . Then 
s 


z(r) 

z(r) 

i x(r) 

x'(r) 

x(r) 

|x' (r) 

^ 0 as  r 


by  virtue  of  (2.3.4).  This  proves  type  (z)  = 0 with  respect  to  X'  . 


Let  z t S and  suppose  type  (z)  = <j  with  respect  to  X . Thus 
lit!  sup  lz(r)/x(r)|  > 0 for  every  x c X . Choose  x'  ^ X'  . Let 

>oo 

s = type  (x')  and  choose  x 6 X^  . Then  the  sequence 


z(r) 

x'(r) 


I z(r) 
|x(r) 


x(r) 

x'(r)| 


satisfies  lim  sup  |z(r)/x'(r)|  >0  by  virtue  of  (2.3.3).  This  proves 

type  (z)  = 0"  with  respect  to  X'  . 

Finally,  let  z C S and  suppose  type  (z)  = s , 0 < s < ct  . Thus 

z(r)/x(r)  0 for  every  x $ X ^ IJ*‘*UX  and  lim  sup  |z(r)/x(r)|  > 0 

s ' 1 CT 

>*oo 

for  every  x € X^  . It  can  be  verified  that  z(r)/x'(r)  ->■  0 for  every 
x'  € X'  ,,U--‘UX'  and  that  lim  sup  lz(r)/x’(r)|  > 0 for  every 

x'  ' using  arguments  similar  to  the  ones  used  in  the  preceding  two 

paragraphs . j^‘| 


CHAPTER  3.  THE  GENERAL  ALGORITHM 


3 . 0 Preliminaries  and  Overview 

In  this  chapter  we  fix  our  attention  on  a particular  linear  differ- 
ence equation 

(3.0.1)  Dy(r)  = g(r)  , r = 0,1,2... 

assuming  that  D t P is  separable  and  g $ 5 . More  specifically, 

Z 

are  interested  in  computing  an  approximation  y of  a particular  solution 

n 

y of  (3.0.1)  which  is  valid  over  some  finite  subsequence  y^’^  . i and 
m given. 

Let  CT  be  the  number  of  distinct  types  in  any  optimally  ranked  basis 
for  D . We  assume  that  y is  known  to  be  a sequence  of  type  t with 
respect  to  D , where  0 5 t < a . We  also  assume  that  initial  values 
y(i)  syCi"*"!) , • • • jyCi'fj”!)  are  known,  where  j is  the  dimension  of  the 
subdominant  subspace  of  type  t for  D . The  extension  of  Theorem  1.3.1, 
from  totally  separable  to  separable  operators,  is  given  in  §3.1.  Accord- 
ingly, these  initial  values  (in  the  presence  of  one  additional  condition) 
suffice  to  determine  the  solution  uniquely. 

Let  k = -^  - j . The  approximating  sequence  y^  is  defined  as  the 
solution  of  (3.0.1)  which  satisfies  the  conditions 

(3.0.2)  y^(i+r)  = y(i+r)  , r = 0,1, . . . , j-1 

and 

(3.0.3)  y (n+r)  = 0 , r = 0,l,...,k-l  . 

n 

r->  • T . T ^ x.Tn 

Since  y is  to  be  approximated  by  y 


, it  is  clear  that  the  value 


of  n tnusc  exceed  m . The  existence  and  uniqueness  of  the  approximating 

sequences  y , for  all  sufficiently  large  n , and  the  convergence  of 
n 

the  sequence  of  aporoximate  values  v (r)  to  y(r)  for  each  value  of  r 

n 

in  the  range  i + j < r < ” , are  proved  under  appropriate  conditions 
in  §3.2. 

The  general  algorithm,  which  is  valid  for  any  separable  operator,  re- 
duces to  the  algorithm  presented  in  §1.3  for  totally  separable  operators. 

We  again  assume  that  the  separable  operator  is  (i,j)  - factorizable,  just 
as  we  did  for  totally  separable  operators.  The  linear  system  (1.3.5)  is 
solved  by  the  factorization  method  by  means  of  a forward  elimination 
stage  followed  by  a back  substitution  stage;  see  §1.3.  The  expansions 
(1.3.19)  and  (1.3.20),  which  are  used  to  estimate  the  optimal  value  of 
n , are  unchanged  for  separable  operators.  In  §3.3  we  discuss  the  stability 
of  the  forward  elimination  stage,  which  was  omitted  in  §1.3.  We  also 
extend  Theorem  1.3.2  to  separable  operators,  showing  that  the  back  substitu- 
tion is  stable.  Finally,  we  note  that  our  remarks  on  the  stability  of 
the  forward  recurrence  solution  of  the  adjoint  equation,  whose  solution 
enters  in  the  truncation  error  expansion  of  §1.3,  apply  to  separable 
operators  as  well  as  totally  separable  operators. 

We  conclude  this  introductory  section  by  fixing  some  notation  for 
later  use  in  this  chapter.  Let 

a 

(3.0.4)  X = U X 

s=l  ® 

be  an  optimally  ranked  basis  for  D , where 


(3.0.5) 


X ={x  ,,x  ,..,,x  },  s = l,2,...,o-. 

s s , 1 s , 2 s , n 

s 


Then  the  dimension  of  the  subdorainant  subspace  of  type  t for  D is 


t 

(3.0.6)  j = I n 

s=l  ^ 


Define  the  set 


(3.0.7) 


t 

U = {u  ,u  , . . . ,u . } = U X 
12  1 s=l  " 


by  means  of  the  correspondence 


X 

p,q 


When  t = 0 it  is  understood  that  j = 0 and  U is  empty, 
define  the  set 


(3.0.8) 


9 ’ ' ■ 


U X : 
s=t+l  ^ 


Similarly , 


when  t = a it  is  understood  that  V is  empty. 


3 . 1 Existence  and  Uniqueness 


Let  y be  a solution  of  (3.0.1)  of  type  t such  that 


(3.1.1) 


y(i+r)  = > r = 0 . 1 , . . . , j -1 


for  some  i > 0 , where  j is  given  by  (3.0.6).  Using  the  notation  |u|( 
for  the  Casoratian  of  U at  i , we  also  suppose  that 


(3.1.2) 


|U| (i)  ^ 0 ; 


see  (3.0.7).  Note  that  for  every  i (3.1.2)  is  satisfied  vacuously  if 
t = 0 , and  also  (by  virtue  of  Casorati's  theorem  [14, §12. 11])  if  t = c 
and  D is  nonsingular.  The  following  theorem  generalizes  Theorem  1.3.1. 
Theorem  3.1.1  If  (3.1.2)  is  valid  then  there  exists  a unique  solution  of 
(3.0.1)  of  type  t such  that  (3.1.1)  is  satisfied. 

Proof:  Suppose  t = 0 . Then  j = 0 and  we  must  prove  that  without  any 

specified  initial  values  y is  uniquely  determined.  Let  v 6 S be  a 
solution  of  (3.0.1),  other  than  y , such  that  type  (v)  = 0 . Since 
V - y € f^(D)  , there  exist  scalars  a $ F such  that 

^ p>q 

n 

CT  p 

y-v=y  Va  X ; 

pil  q=l  P’P 


see  (3.0.5).  At  least  one  of  the  a is  nonzero;  let  the  largest 

p.q 

value  of  p such  that  a 0 for  some  q,  l5q5n,be  p = s. 

P.q  p 

Then 


lim  y (r) -V ( r) 

r-rco  X (r) 
s , 1 


because  type  (y)  = type  (v)  = 0 . By  part  (i)  of  Definition  2.2.1,  we 
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have 


^ . s-1  n a X (r) 

lim  r P,q  P»q 


p=l  q=l  s,l 


= 0 


It  follows  from  the  above  three  equations  that 


n 


lim 

>oo 


s a X (r) 

I = 0 

X,  1 (r) 

q-1  s , 1 


But  this  result  contradicts  part  (iii)  of  Definition  2.2.1;  therefore  we 
conclude  that  y = v . 

Alternatively  suppose  t > 0 . Let  v ^ S be  a sequence  of  type  t 
which  satisfies  (3.0.1)  and  is  coincident  with  y(i+r)  for  r = 0,1,..., 
j-1  . Then  there  exist  scalars  a 6 F such  that 

p>q 

n 

V = y + I I a „ • 

p=l  q=i  P'1  P'<! 


Let  the  largest  value  of  p in  this  sum  such  that  at  least  one  of  the 
scalars  a is  nonzero  be  p = s . Then  s ^ t , because  s > t 

p,q 

would  imply 


lim  v(r)-y (r) 

r-«°  x , (r) 

S j 1 


0 


(since  type  (v)  = type  (y)  = t)  and 


lim 

>oo 


S*1 

I 


p=l 


pa  X (r) 

y p>q^>q 


0 


(by  part  (i)  of  Definition  2.2.1).  This  gives 


lim 


= 0 , 


s a X (r) 
y s.q  s,q 

q-1  s,l 


which  is  a contradiction. 

Thus,  a =0  whenever  p > s . The  a having  p 5 s satisfy 

p»q  p,q 


the  system  of  linear  equations 


s n 


y y a X (i+r)  = 0 , r = 0,1,..., j-1  . 

pil  qii  p.q  p,q 


This  system  is  nonsingular  by  the  condition  ju|(i)  7^  0 , so  we  have 

G =0  for  all  p and  q . Thus  v(r)  = y(r)  for  all  r and  the 

p,q 

theorem  is  proved.  [TT] 

Theorem  3.1.2  If  D is  nonsingular  then,  for  each  i > 0 , !u|(i+r)  ^ 0 

for  at  least  one  value  of  r in  the  range  0 5 r 5 -2  - j . 

Proof:  Since  D is  nonsingular,  Casorati's  theorem  is  valid,  i.e., 

iX'(i)  f 0 for  all  i . Theorem  3.1.2  is  then  an  immediate  consequence 
of  Laplace's  general  theorem  on  the  expansion  of  determinants  [15,  §93].  Q 
Theorem  3.1.2  implies  there  is  no  shortage  of  points  i which  satisfy 
condition  (3.1.2)  when  the  difference  operator  is  nonsingular.  Our  next 
goal  is  to  show  that  any  optimally  ranked  basis  may  be  used  to  locate  an 
admissible  value  of  i . 

Theorem  3.1.3  If  |uj(i)  r 0 for  some  fixed  i > 0 , then  the  leading 
principal  minor  of  order  j of  the  Casoratian  evaluated  at  i of  every 
optimally  ranked  basis  for  D is  nonzero. 

Proof:  If  j = 0 then  there  is  nothing  to  prove,  so  let  us  assume  j > 0 . 

Also,  let  us  assume  temporarily  that  ] < I . Let  X'  be  an  optimally 
ranked  basis  distinct  from  X , where 


X’ 


o 

U X’ 

1 s 
s=l 
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and 


XL  = {x'  ,x’  , 

S s, 1 s,2 


,x  j-  , 
s ,n 

s 


s 1,2,.. .jO”  \ 


compare  (3.0.4),  (3.0.5).  Define  the  sets 


U'  = {u' ,u ' . . . ,u' } = U X' 
^ 3=1 


and 


V = {v’  v’  = u X'  . 

^ ^ s=t+l  = 


There  exist  scalars  a € F such  that 

p.q 


P qil  P>P  P ,=5+1  P-'l  P-2 


and 


j 

v'.=  ya  u + / a V. 

P-3  qii  P,q  q P,q  q-3 


p = j-l-1,  j+2, . . . ,^  . 


The  matrix  of  coefficients  (a  ) is  nonsingular  because  both  X and  X' 

P»  q 

are  bases  of  the  finite-dimensional  linear  space  K (D)  . Furthermore, 


from  the  definition  of  an  optimally  ranked  basis  it  is  easy  to  show  that 


(a  ) is  block  lower  triangular  with  the  s-th  diagonal  block  having  order 

p>q 

n . Because  of  the  block  triangular  structure,  the  determinant  of  (a  ) 

s p,q 

is  equal  to  the  product  of  the  determinants  of  the  diagonal  blocks.  And 


because  the  determinant  of  (a  ) is  nonzero,  so  is  the  determinant  of 

P.q 

every  diagonal  block.  Therefore,  the  leading  principal  minor  of  order  j 


of  the  determinant  of  (a  ) is  nonzero.  Let  A.  denote  this  minor. 

P.q  3 


Next,  consider  the  Casorati  matrix  [X](i).  It  can  be  written  in 


block  form  as 


[X](i)  = 


[U](i) 


[V]  (i+j) 


where  the  actual  form  of  the  blocks  indicated  by  asterisks  is  unimportant 
Similarly,  we  write 


[X'](i)  = 


[U'](i) 


[V'](i+j) 


If 


(a  ) = 
P i 1 


is  the  corresponding  block  form  of  (a  ) , we  verify  readily  that 

P ,q 


A 

0 

11 

A , 

A „ 

21 

22 

[U](i) 


[X'](i)  = 

[V](i+j) 

The  leading  principal  minor  of  order  j of  |X'j(i)  is  evidently  given  by 


^11 

•A^n 

. 21 

0 

• Tr 

:^22 

U'|(i)  = |U|(i) -det(A^p  . 


Since  |Uj(i)  ^ 0 , by  hypothesis,  and  det(A^^)  = 0 we  have 

!U'|(i)  r 0 . This  completes  the  proof  for  j < -3  . 

Finally,  we  note  that  the  proof  in  the  case  j = ^ follows  by  an 
argument  analogous  to  Chat  just  used  for  j < £ , except  that  there  is 
no  partitioning  of  matrices  since  V and  V are  empty.  j'_J 
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3 . 2 Convergenc  e 

y 

Assume 

(3.2.1)  iu|  (i)  0 

and  let  y be  a solution  of  type  t of  (3.0.1),  where  0 5 t 5 ct  . For 

each  n > m + 1 such  that  a solution  y of  (3.0.1)  exists  satisfying 

n 

both  (3.0.2)  and  (3.0.3),  there  exist  scalars  (n)  ,U2(n)  , . . . ,a^  (n)  and 
j3^(n)  ,p^(n)  , . . . ,p^(n)  such  that 

j k 

(3.2.2)  ^ ^ ^ ^ ' 

S=1  S=1 

In  the  case  t = a , we  have 

y = y+  V a(n)u 
^-i  s s 

s=l 

Then  we  prove  readily,  using  (3.0.2)  and  (3.2.1),  that  ~ ^ 

s = 1,2,...,^  . Thus  for  each  n , v (r)  = v(r)  for  r = i.i+1 

n " 

Therefore  we  restrict  the  ensuing  discussion  to  the  cases  in  which 
0 5 t < c . 

Consider  the  set  C of  all  subsets  of  k distinct  sequences  from 

X . We  have  k > 0 by  our  assumption  that  t < c . Obviously, 

V = {v  ,v  , . . . ,v  } € C . We  shall  say  that  the  optimally  ranked  basis  X 
12  k 

is  j -normal  provided  that  |v| (r)  ^ 0 for  all  sufficiently  large  r , and 

(3.2.3)  |v|(r)  =o{|v|(r)}  as  r ->  «> 

for  every  V 6 C - {V}  . When  t = 0 , C = {V}  and  X is  automatically 


^The  reader  is  referred  to  §3.0  for  notation  and  underlying  assumptions. 


j-normal  (0-normal) , When  t > 0 the  sets  V contain  at  least  one 


sequence  from  U = {u^ , u^, . . . , u^ } . In  this  case  (3.2.3)  expresses  the 
quite  reasonable  condition  that  whenever  one  or  more  of  the  "dominant" 
solutions  V, ,v^,...,v,  is  replaced  in  the  Casoratian  |vl  by  a "sub- 
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dominant"  solution,  then  the  resulting  Casoratian  is  dominated  by  |v|(r) 
as  r . 


We  now  state  and  prove  the  convergence  theorem.  A similar  result, 
but  only  for  homogeneous  linear  difference  equations,  is  given  by  Zahar 
[27 , Th.  5.1] . 

Theorem  3-2.1  Assume  X is  a j-normal  optimally  ranked  basis  such  that 
(3.2.1)  is  satisfied  and  y(r)  is  a solution  of  (3.0.1)  of  type  t.  Also 


assume 


(3.2.4) 


1 (y) I (r) 

lim  ' s 


Vi  (r) 


=0,  s=l,2,...,k 


where  V (y)  = {v  , .,,,v  ,y,v  ,...,v  } . Then  for  sufficiently  large 

S J-  S ' i. 

n there  exists  a unique  solution  y of  (3.0.1)  which  satisfies  (3.0.2) 

n 

and  (3.0.3).  Furthermore, 


(3.2.5) 


lim 

n-v<» 


y (r) 

n 


y(r) 


for  each  fixed  value  of  r > i + j . 

Proof:  From  (3.2.2),  (3.0.2)  and  (3.0.3)  we  derive  the  linear  system 


U^(i)  ... 

(i) 

v^(i)  ... 

V^(l) 

a^(n) 

“ 0 

u^(i+j-l).  . 

. u . (i+j-1) 

(i+j-1).. 

. Vj^(i+j-L) 

.bl"?. 

0 

Uj^(.n) 

u . (n) 

J _ 

v^(n)  ... 

v^(n) 

P^(n) 

-y(n) 

u ( n+k- 1 ) . . 

. u^ (n+k-1) 

v^ (n+k-l) . 

.Vj^(  n+k-1) 

c 

CO. 

1 

-y ( n+k-1) 

Although  the  block  structure  indicated  in  this  system  disappears  when 
t = 0 , this  modification^ causes  no  difficulty  in  the  following  proof. 

Suppose  V were  such  that  the  upper  right  block  of  the  matrix  is 
zero  (again,  this  is  satisfied  vacuously  if  t = 0) . Then  the  linear 
system  becomes 


U^(i)  ... 

(i) 

0 

0 

a^(n) 

0 

u^(i+j-l) . . 

. u . (i+j-1) 

0 

0 

a.  (n) 

0 

u^(n) 

Uj  (n) 

v^(n) 

Vk(n) 

P^(n) 

-y(n) 

u^(n+k-l) . 

.u . (n+k-1) 
J 

v^  (n-hk-1)  . 

. .V  (n-hk-1) 
k J 

-y  ( n-hk- 1 ) 

The  determinant  of  this  system  is  | U j (i) * | V | (n)  , and  it  is  nonzero  for 
all  sufficiently  large  n , in  view  of  our  assumptions.  Using  Cramer's 
rule  and  Laplace's  general  theorem  on  the  expansion  of  determinants 
[15, §93],  we  then  find  that 


^iCn)  = a^(n)  = • • • = (n)  = 0 

and 


P^(n) 


Vg(y)  I (n) 
jvj  (n) 


s = 1,-2, 


,k  . 


Therefore,  for  all  r and  for  all  sufficiently  large  n , we  have 


(3.2.6) 


y^(r)  = y(r) 


k 

r 

L 

s=l 


I V^(y)  I (n) 


V (r)  . 


The  other  case  in  which  the  block  structure  disappears,  which  is 
t = c , was  disposed  of  earlier. 


Then  (3.2.5)  follows  from  (3.2.4). 


In  order  to  complete  the  proof,  let  us 

(3.2.7)  v'  = V + r u,  + + Y .u 

P P P,1  1 P,J  j 


introduce  the  sequences 
, p = 1,2, ... ,k  , 


where  the  y are  to  be  chosen  in  such  a way  that 

p.q  ^ 

v^(i+r)  =0  , r = 0,1,..., j-1  . 

Thus  the  y have  to  satisfy  the  linear  systems 

p>q 


u (i) 

u.(i) 

y , 

-V  (i) 

1 

u^(i+j-l) . . 

J _ 

.u. (i+j-1) 

P,1 

= 

P 

-V  (i+j-1) 

J J 

^ P’jJ 

L P J 

each  of  which  is  nonsingular  because  |u|(i)  0 , by  hypothesis.  Let 

V'  = {v' , v' , . . . , v' } . Clearly  it  suffices  to  show  that  U U V is  a 
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j -normal  optimally  ranked  basis  for  D such  that 


(3.2.8) 


r->-«  I V I ( r ) 


= 0 


s = 1,2, 


,k  . 


Since  v dominates  every  sequence  in  U , we  have  from  (3.2.7) 
P 


v’(r) 

^ . -^1  as  r ->-  «> 
V (r) 

P 


for  each  p , 1 < p < k . Therefore  v'  has  the  same  type  as  v 

P P 

Furthermore,  we  may  verify  the  three  conditions  of  Definition  2.2.1, 
as  follows: 


i)  ■•‘(f)  ^ u(r)  ^ 

v'(r)  v(r)  v'(r) 
P P P 

for  every  u € U , and 


0 as  r 


oo 
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v'(r)  v'(r)  V (r)  v (r) 

— 2_ = —2— . — 2_ . _2_ ->Q  og  1- 

v'(r)  v(r)  v(r)  v'(r) 

a p q q 


whenever  type  (v')  < type  (v'); 

P q 


ii) 


v'(r)  v'(r)  V (r) 

_2 2 . _2 


V ' (r)  V (r)  V (r) 

q p q 


v'(r)  V (r) 

q q 


as  r “ 


whenever  type  (v^)  = type  (v'); 

(iii)  Let  v',v'  ,...,v'  6 V be  all  the  sequences  of  one  type, 

P P+1  q 

and  let  v'  be  any  fixed  one  of  these.  Let  the  corresponding  unprimed 

sequences  be  v,v  ,v  v 6 V . Take  5 ,5  S € F , not  all 

P P+1  q P P+1  q 

zero.  For  each  s = p,p+l,...,q  we  have 


5 v'  5 (v  +Y  ,u  +---+V  .u.)  6 V 

s s _ s s s,  1 1 3,1  1 _ s s ... 


V 


(since  are  of  lower  type  than  v’  = v{l+o(l)})  . Hence 


Furthermore 


6v’+***+5  v'  S V +---+6  V 

= . 

V V 


5 v'+*’’+6  v'  S V +‘*-+6  V 


P .P 


- P.  P- 


q q 


V 


{1+0(1)}  + 0(1) 


from  which  it  is  clear  that 


lim  sup 

■J-'— woo 


I 6 v'(r)/v'(r) 
s s 

s=p 


> 0 


Thus,  U U V is  an  optimally  ranked  basis. 


Next  we  prove  that  |v'j(r)  ^ 0 for  all  sufficiently  large 


r . 


We  recall  that  a determinant  is  a multilinear  function  of  its  columns. 
If  we  express 


V'(r) 


IV 


1’ 


’ V 

2’ * ’ •’  k 


(r) 


where  v'  v'  ...,v'  temporarily  denote  the  successive  columns  of  | V | (r)  , 
12  ic 

then  we  derive 


|V'|(r)  = lv^,v^,...,v^|(r)  + 'f'l  ,q  i ""q ‘ 


q=l 


= V 


,v  ,v'  . . . ,vj  (r)  + I r ju  ,v'  . . . ,vj  (r) 


1’  2’  3 


a=l 


l,q'  q’  2- 


and  hence 


|v'^(r)  = |v!(r)  + I T,  | u , v ' . . . , v ' | ( r) 


q=l 


l,q'  q’  2 


q=l 


+ • • • 


q=l 


h.q'''r''2 


If  we  continue  this  construction  until  every  reference  to  a primed  v is 
removed,  then  it  is  clear  that  |v'|(r)  is  equal  to  the  sum  of  |v|(r) 
plus  a linear  combination  of  determinants  formed  from  | V j ( r)  by  replace- 


ment  of  at  least  one  v by  a u . By  our  assumption  that  X is  j -normal, 
we  conclude 

(3.2.9)  iV'l(r)  = |vj  (r){l+o(l)}  as  r . 

Since  !vj(r)  ^ 0 for  sufficiently  larger,  we  also  have  |V  | (r)  ^ 0 
for  sufficiently  large  r. 

If  one  or  more  of  the  solutions  replaced  in 

I V j (r)  by  subdominant  solutions  from  U , then  a construction  like  the 
one  preceding  shows  chat  the  modified  Casoratian  |v'j(r)  is  dominated  by 
|v|(r)  as  r ® , since  the  original  optimally  ranked  basis  is  j-normal. 
Consequently,  using  (3.2.9),  we  conclude  U € V is  j-normal.  A similar 
argument  shows  that  (3.2.8)  is  satisfied,  and  the  theorem  is  proved, 


3.3  Stability 


Let  us  consider  solutions  of  (3.0.1)  of  type  t,  where  0 5 t < ex  . ' 

Let  j be  the  dimension  of  the  subdominant  subspace  of  type  t for  D . 

Let  i be  a point  at  which  |u!(i)  0 ; see  (3.0.7).  Suppose  D is 

(i, j)-factorizable;  see  Definition  1.2.1.  Let  = [a^, a^, . . .a^]  be  a 

lert  factor  and  B ^ = [b^  ^ , . . . , b^^  '^  ] the  corresponding  right 

factor  of  , so  that  ^ is  an  (i, i)-factorization  of  D . 

3 ' 3 3 

The  next  theorem  provides  a generalization  of  Theorem  1.3.2,  from 
totally  separable  linear  difference  equations  to  equations  that  are  merely 
separable.  The  remark  following  Theorem.  1.3.2  applies  here  as  well.  That 
is,  the  back  substitution  stage  of  the  algorithm  is  stable,  at  least  for 
sufficiently  large  r. 

Theorem  3.3.1  If  lu|(i)  7^  0 and  = A^B^"*"^  , then  the  difference 
operator  is  separable  and  any  solution  y of  (3.0.1)  of  type  t or 

less  is  a sequence  of  type  zero  with  respect  to  . 

Proof:  Suppose  3=0,  i.e.,  type  (y)  = 0 with  respect  to  D.  Then 

A^  = [a^]  and  B^  = [b^ , b^ , . . . , b^]  and  we  have  = A^B^  . 

0 0 1 V 

Now  let  X € X , where  X is  the  optimally  ranked  basis  (3.0.4).  Then 
Dx  = 0 ; hence  D^x^  = A^B^x^  = 0 . This  last  equation  is  equivalent  to 

ap,(r)  y b (r)x(r+s)  = 0 . r = i,i+l,...  . 

0 ■'n  s 

s=0 

Since  a^  is  free  of  zeros  by  the  definition  of  (i,j)  - factorizability , 
we  have  = 0 . Thus  {x^|x  i X}  is  an  optimally  ranked  basis  for 

and  it  follows  im.mediately  that  type  (y)  = 0 with  respect  to  B 

The  case  t = o corresponds  to  a maximal  solution  of  (3.0.1).  Since  th 
general  algorithm  reduces  to  pure  forward  recurrence  in  this  case,  it  is 
stable ; see  ^1  . 1 . 
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Next  suppose  j satisfies  0 < j < -S  . Let  v 6 V , where  V is 
given  by  (3.0.8).  In  view  of  our  assumption  that  ju|(i)  0 , we  may 
assume  that 


v(i+r)  =0,  r=0,l,...,j-l; 

compare  the  proof  of  Theorem  3.2.1.  Since  v £ /((D)  we  have  Dv  = 0 . 

It  follows  that  for  each  n > i + j + max(j,k)  the  subsequence  ^ 

satisfies  the  finite  boundary  value  problem 


pi,n-j-l^i+j ,n-l 
3 


p 

^n-'^  ,n-j-l^n,  n+k-1 


compare  (1.2.4).  For  each  n,  these  boundary  value  problems  are  nonsingular 
and  factorizable  as  shox^m  by  Theorem  1,2.1. 

Therefore,  Theorem  1.2.2  is  applicable  and  we  see  that  ^ satisfies 

i+j,n-l  i+j,n-l  _ 

»0 

for  each  n > i + j + max(j,k)  ; compare  (1.2.7)  and  the  remark  following 
the  proof  of  Theorem  1.2.2.  The  last  equation  is  equivalent  to 

k 

I b (r)v(r+s)  = 0 , r = i+j . i+j+l , . . . .n-1  . 
s=0  ® 


0 

„n-k,n-l  n,n+k-l 
B,  V 


Since  n may  be  arbitarily  large,  we  have 
k 

I b (r)v(r+s)  = 0 , r = i+j , i+j+1 , . . . 
s=0  ^ 

or,  equivalently,  = 0 . Therefore,  {v^"*"^|v  6 V}  is  an  optimally 

ranked  basis  for  B^"^^  . Since  y is  dominated  by  every  sequence  in  V , 


ic  follows  that  type  (y)  = 0 with  respect  to  . | | 

Now  let  us  turn  to  the  for^^7ard  elimination  stage  of  the  algorithm. 
This  consists  of  applying  a finite  number  of  steps  of  Gaussian  elimina- 
tion, without  pivoting,  to  the  infinite  linear  system 


(3.3.1) 


i i+i  i 

D.y  ^ = g - 


y 


compare  (1.3.6).  Rounding  errors  introduced  in  this  process  are  equivalent 
to  perturbations  in  the  original  problem.  Thus,  instead  of  satisfying 
(3.3.1)  the  computed  solution  y^  ^ , say,  is  an  exact  solution  of  a 
system  of  the  form 


(3.3.2) 


(D^  + = g"*  - 


•^i,i+j-y,i+j-i 

b 


+ e' 


where  the  terms  and  e^  represent  the  perturbation. 

The  effect  of  introducing  a single  element  of  e(r)  , say  at  r = s , 
is  to  perturb  the  true  solution  y^"^^  by  a linear  combination  of  the 
solutions  of  the  corresponding  homogeneous  equation.  Because  of  the 
boundary  conditions  we  have  imposed  "at  infinity",  only  the  solutions 
u^.u^j.-.u^  enter  from  the  optimally  ranked  basis  {u^.u^ , • • . , u^  , v^  , v^ , 
...,v,  } ; compare  equation  (3.0.3)  and  Theorem  3.2.1.  Subsequently, 

eC 

we  have  only  to  consider  what  happens  when  r > s , and  here  it  will  be 

the  multiple  of  u^  that  will  be  the  fastest  growing. 

The  effect  of  , also,  can  be  allowed  for  by  making  a perturbation 

3 

of  the  right  side.  To  terras  of  the  first  order,  we  have 


73 


Then  the  arguments  of  the  preceding  paragraph  again  apply. 

Since  the  wanted  solution  y is  of  type  t , we  conclude  that  the 
process  of  forward  elimination  is  stable  in  the  sense  that  each  perturba- 
tion subsequently  grows  at  a rate  that  does  not  exceed  that  of  che  wanted 
solution.  Of  course,  there  will  be  a loss  of  accuracy  when  an  element  of 
(or  e^)  is  large  compared  with  the  rounding  error  in  the  stored  value 
of  the  corresponding  element  of  (or  of  the  right  side  of  (3.3.1)). 

This  may  occur,  for  example,  when  there  is  heavy  cancellation  in  the 
formation  of  an  element  that  is  subsequently  used  as  a pivot.  Hoi>;ever, 
the  important  point  is  that  the  effect  of  each  such  loss  is  not  magnified 
in  subsequent  steps  of  the  algorithm. 

An  extension  of  the  foregoing  discussion  shows  that  if  we  select  a 
value  of  j (that  is,  the  number  of  prescribed  initial  conditions)  that  is 
less  than  the  dimension  of  the  subdominant  subspace  of  the  wanted  solution, 
then  the  forward  elimination  remains  stable.  However,  the  back-substitution 
is  now  unstable;  compare  the  proof  of  Theorem  3-3.1.  Similarly  if  j is 
too  large  then  the  forward  elimination  is  unstable,  and  the  back-substitu- 
tion is  stable.  Nevertheless  if,  in  practice,  the  actual  instabilities  are 
weak  owing  to  weak  separation  of  solutions  of  consecutive  types  as  r , 

then  it  may  be  advantageous  to  select  a value  of  j that  differs  from 
the  dimension  of  the  subdominant  subspace  of  the  wanted  solution.  This  is 
because  the  loss  of  accuracy  caused  by  the  instability  is  offset  by  a 
substantial  reduction  in  the  value  of  the  terminal  point  n owing  to 
increased  convergence.  This  modification  is  illustrated  by  numerical 
examples  in  Chapter  4. 


CHAPTER  4.  NUMERICAL  EXAMPLES 


i . 0 Introduction 

Several  numerical  examples  will  be  given  to  illustrate  the  general 
algorithm.  The  examples  all  involve  fourth-order  linear  difference 
operators  with  nonconstant  coefficients.  Both  homogeneous  and  inhomo- 
geneous difference  equations  are  treated. 

The  fourth-order  operators  were  produced  by  using  the  method  described 
in  [14,  §12.22],  starting  with  two  second-order  recurrence  relations  which 
have  as  solutions  Bessel  functions,  modified  Bessel  functions  or  associated 
Legendre  functions.  The  inhomogeneous  equations  have  as  particular 
solutions  Anger-Weber  functions  or  Struve  functions.  The  right  sides  of 
these  equations  were  produced  as  a by-product  of  the  method  used  to  pro- 
duce the  fourth-order  operators. 

The  production  of  the  fourth-order  difference  operators  from  the 
corresponding  pairs  of  second-order  operators  requires  only  elementary 
algebraic  manipulations.  These  calculations  are  rather  lengthy,  however, 
and  would  be  extremely  difficult  (as  well  as  tedious)  to  complete  accurately 
by  hand.  Instead,  the  MACSYMA  symbol  manipulation  system,  described  in 
[12],  was  employed  . This  procedure  had  the  added  advantage  of  storing 
the  required  formulas  for  the  coefficients  directly  in  the  computer,  ready 
for  numerical  evaluation,  thereby  further  reducing  the  possibility  of 
human  error.  Similarly,  MACSYMA  was  used  to  advantage  in  the  production 
of  the  right  sides  for  the  inhomogeneous  examples. 


Development  of  which  is  currently  supported,  in  part,  by  the  United  States 
Department  of  Energy  under  Contract  Number  DE-AC02-79ER10357  and  by  the 
National  Aeronautics  and  Space  Administration  under  Contract  Number  NSG1321. 


In  the  case  of  the  Bessel  functions  J (x)  , Y (x)  and  modified 

r r 

Bessel  functions  I (x)  , K (x)  we  restrict  ourselves  to  integer  order 

r r 

r > 0 and  real  argument  x > 0 . Relevant  properties  of  these  functions 
are  given  in  [16,  Chapter  9].  For  example,  for  fixed  x the  functions 
J (x)  and  Y (x)  satisfy  the  linear  recurrence  relation 

r - ^ ' 


(4.0.1) 


y(r-l)  - — y(r)  + y(r+l)  = 0 


Similarly,  I^(x)  and  (-)  K^(x)  satisfy 


(4.0.2) 


y(r-l)  - — y(r)  - y(r+l)  = 0 


The  Anger -Weber  functions  E^(x)  satisfy 


(4.0.3) 


y(r-l)  - ^ y(r)  + y(r+l)  = 


2{l-(-l)  } 


and  the  Struve  functions  H (x)  satisfy 

r 


(4.0.4) 


y(r-l)  - ^ y(r)  + y(r+l)  = 


,1  . r 

x) 

VnF(r+  |) 


compare  (4.0.1)  in  both  cases..  Relevant  properties  of  the  E^(x)  and 
H^(x)  appear  in  Giver's  paper  [21]  where  they  were  used  as  examples  for 
second-order  inhomogeneous  equations. 

The  coefficients  of  the  fourth-order  homogeneous  recurrence  relation 
obtained  from  (4.0.1)  with  x = x^  and  (4.0.2)  with  x = x^  are  shown 
in  Figure  6 in  the  form  of  FORTRAN  statements  produced  by  MACSYMA.  In 


addition.  Figure  6 shows  the  right-hand  sides  corresponding  to  E^(x^) 
and  H^(x^)  . We  designate  the  difference  operator  D defined  by  these 
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coefficients  the  JYIK  operator . Since  the  functions  ’ 

^r^^2^  , (-)  K^Cx^)  are  linearly  independent,  they  form  a fundamental 
set  of  solutions  of  Dy  = 0 , 

Similarly  in  the  case  of  the  associated  Legendre  functions  P’j^(x)  , 

Q^(x)  we  restrict  consideration  to  integer  degree  r > 0 , integer  order 
u > 0 , and  real  argument  x > 0 . Relevant  properties  of  these  functions 
may  be  found  in  [16,  Chapter  8].  We  use  the  r-wise  linear  recurrence 
relation 

(4.0.5)  (r+u)y(r-l)  - (2r+l)xy(r)  + (r-u+l)y (r+1)  = 0 , 

valid  for  fixed  x and  u . The  coefficients  of  the  fourth-order  differ- 
ence operator  D,  shown  in  Figure  7,  were  produced  (using  MACSYMA)  from 
(4.0.1)  with  X = x^  and  (4.0.5)  with  x = x^  . Figure  7 also  shows  the 
right  sides  corresponding  to  E^(x^)  and  H^(x^)  . In  this  case  the 

functions  J (x^ ) , Y (x^ ) , p'^(x.)  , Q^(x„)  form  a fundamental  set  of 
r 1 r 1 r 2 r 2 

solutions  of  Dy  = 0 . Accordingly,  we  designate  D the  JYPQ  operator. 

A FORTRAN  program  has  been  written  to  implement  the  general  algorithm. 
The  user  may  select  either  the  JYIK  operator  or  the  JYPQ  operator.  He  may 
also  elect  to  solve  the  homogeneous  equation,  the  inhomogeneous  equation 
having  the  Anger-Weber  function  as  a solution,  or  the  inhomogeneous 
equation  having  the  Struve  function  as  a solution.  Alternatively,  he  may 
provide  a subroutine  which  generates  the  coefficients  of  an  arbitrary 
linear  difference  equation  of  any  order. 

The  user  supplies  the  program  with  the  starting  point  i and  the  terminal 
point  m of  the  desired  subsequence  of  the  solution  to  be  computed,  the 
number  j of  initial  values,  the  subsequence  y^’^***^  ^ of  initial  values, 
and  the  parameters  v and  e for  use  in  the  convergence  test  (to  be  described 
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below) . The  program  solves  the  proposed  problem  by  executing  three  phases 

In  the  initial  phase  the  first  m+l+v-i-j  rows  of  the  infinite 

linear  system  (1.3.6)  are  stored.  Forward  elimination  without  pivoting 

is  performed  on  these  rows.  This  results  in  the  storing  of  the  first 

m+l+v-i-j  rows  of  the  left  and  right  factors  of  the  (i,j)~ 

factorization  ^ , and  also  the  subsequence 

J J 

solution  of  (1.3.8).  The  elimination  procedure  is  carried  out  by  modifica 
tions,  described  below,  of  the  subroutines  DGBFA  and  DGBSL  which  are 
available  in  the  mathematical  software  package  known  as  LINPACK;  see  [ 8 ] . 

The  unmodified  DGBFA  performs  forward  elimination  with  partial 
pivoting  on  a finite  band  matrix,  ignoring  the  right  side.  It  produces 
an  upper  triangular  matrix,  a lower  triangular  matrix  and  a vector  of 
pivoting  information.  These  results  may  then  be  used  by  the  unmodified 
DGBSL  to  solve  the  corresponding  linear  system  with  specified  right  side. 
Alternatively,  the  user  can  direct  that  DGBSL  solve  a linear  system  having 
as  its  matrix  of  coefficients  the  transpose  of  the  matrix  input  to  DGBFA. 
Provision  has  been  made  in  DGBFA  to  suppress  the  partial  pivoting.  When 
this  is  done  the  vector  of  pivoting  information  is  superfluous  and  (except 
for  rounding  errors)  the  matrix  obtained  by  premultiplying  the  upper 
triangular  matrix  by  the  lower  triangular  matrix  is  the  original  band 
matrix.  This  factorization  is  therefore  of  exactly  the  kind  required  by 
Definition  1.2.1;  see  equation  (1.2.5). 

When  partial  pivoting  has  been  suppressed  in  DGBFA,  DGBSL  proceeds, 
in  effect,  by  forward  elimination  without  pivoting  followed  by  back 
substitution.  Provision  has  been  made  in  DGBSL  to  exit  after  the  forward 
elimination  stage  and  to  re-enter  at  the  start  of  the  back  substitution 
stage.  Thus  the  user  has  the  option  of  examining  the  result  of  the 
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forward  elimination  stage  and  perhaps  electing  to  continue  it  without 
having  to  do  any  back  substitution.  All  of  the  above  modifications  of 
DGBFA  and  DGBSL,  to  provide  just  the  tools  that  are  needed  in  our 
algorithm,  were  easily  achieved  because  of  the  clarity  and  flexibility  of 
the  coding  and  documentation  of  LINPACK. 

We  return  to  the  description  of  the  initial  phase.  The  next  computa- 
tion is  the  subsequence  . compare  equations  (1.3.16)  and 

m 

(1.3.17)  with  s = m . Since  w^  satisfies  the  adjoint  equation,  it 

m 

is  produced  using  the  forward  elimination  stage  of  the  (modified)  sub- 
routine DGBSL  in  which  the  input  matrix  to  DGBFA  is  automatically  trans- 
posed. Finally,  the  initial  phase  computes 

r-1 

(4.0.6)  y (m)  = V w (s)z(s)/b  (m) 

r ^ m o 

s=m 

for  r = m+1 ,m+2 , . . . ,m+v ,m+v+l  and 

r+v-1 

(4.0.7)  r\  (m)  = 7 w (s)z(s)/b  (m) 

r,v  ^ m o 

s=r 

for  r = m+1  . Equation  (4.0.6)  is  the  same  as  (1.3.18)  apart  from  a 
change  of  symbols.  Equation  (4.0.7)  is  the  approximation  to  the  truncation 
error 

Tj^(m)  = y(m)  - y^(m) 

which  is  obtained  when  y(m)  is  replaced  by  y , (m)  : compare  (1.3.20). 

r+'j 

After  completion  of  the  initial  phase,  the  program  proceeds  to  the 
iteration  phase.  One  additional  row  of  the  infinite  linear  system  (1.3.6) 
is  stored.  Since  the  previously  determined  part  of  the  ( i, j ) -factorization 
has  replaced  the  corresponding  rows  of  the  original  matrix  in  the  computer 


storage,  these  rows  are  no  longer  available  to  extend  the  forward  elimina- 
tion stage.  Therefore,  extension  of  the  forward  elimination  is  accomplish- 
ed by  an  alternative  procedure  that  requires  the  solutions  of  two  triangu- 
lar linear  systems  of  orders  j and  k . This  is  achieved  using  DTRSL,  a 
UNPACK  subroutine  that  solves  triangular  systems.  Next,  the  formerly 

obtained  subsequences  of  the  solutions  and  w™"*"^  are  extended  to 

m 

include  the  next  higher  term.  Then  the  series  (4.0.6)  and  (4.0.7)  are 
summed  with  the  new  value  of  r , and  the  following  convergence  test  is 
applied:  If 


(4.0.8) 


then  the  program  proceeds  to  the  third  and  final  phase;  otherwise,  the 
iteration  phase  is  repeated  and  the  convergence  test  re-applied. 

The  final  phase  of  the  program  is  the  back  substitution  stage  of  the 
algorithm.  Let  n denote  the  value  of  r which  satisfies  (4.0.8).  The 

back  substitution  is  begun  at  this  point,  using  the  (modified)  subroutine 
DGBSL.  This  yields  the  computed  solution  ^ (1.3.7)  and  also 

the  desired  approximation  of  y(r)  over  the  range  r = i+j , i+j+1 , . . . ,m  . 

In  all  of  the  examples  presented  in  this  thesis  the  chosen  termination 

criterion  corresponds  to  an  estimated  relative  error  of  10  in  the  final 

term  y (m)  of  the  requested  sequence  of  values  of  the  function.  Thus 
n 

t = 10  in  equation  (4.0.8).  For  the  values  of  x^  and  x^  selected 
in  the  examples  the  solutions  are  reasonably  well-separated  as  r , 

hence  we  take  the  "testing  parameter"  v in  equation  (4.0.7)  to  be  1 
throughout.  Initial  values  were  entered  to  10  significant  figures  in 


every  case.  All  computations  were  performed  internally  to  18  figures 


and  the  output  printed  to  10  figures.  For  each  example  we  present  (i) 

the  solution  z (i+j  ) , z (i+j-fl)  , . . . ,z  (n)  of  equation  (1.3.8)  (the  forward 

elimination  stage);  (ii)  the  solution  w (m+1)  (m+2)  , . . . ,w  (n)  of 

mm  m 

(1.3.16),  (1.3.17)  (the  forward  recurrence  of  the  adjoint  equation  for  use 

in  the  termination  criterion);  (iii)  the  sequence  y ^^(m),y  .(m),..., 

m+1  m+2 

y (m)  of  approximants  to  the  value  y(m)  (see  equation  (4.0.6)); 
n“K; 

(iv)  the  sequence  ri  (m)  ,"n  , » , (m)  , . . . ,ri  (m)  of  estimates  of  the 

m+1 , 1 m+  2,1  n , 1 

truncation  error  incurred  in  accepting  the  terminal  points  of  m+1 ,m+2, . . . , n 
respectively  (see  equation  (4.0.7));  (v)  the  sequence  y^(n-l) , y^(n-2) , . . . , 

y^(i+j)  derived  from  equation  (1.3.7)  (the  back  substitution  stage). 

The  initial  values  that  were  used  in  the  boundary  value  problem  actu- 
ally solved  in  each  example  are  shox«m  in  parentheses  in  the  tables;  see, 
for  example.  Table  5.  The  terminal  values  (which  are  all  zero)  are  also 
shown.  Forward  recurrence  of  the  j-th  order  difference  equation  A^z^  = g^ 
starting  from  the  given  initial  values  produces  a solution  of  (1.3.8)  (the 
algorithm  produces  this  solution  by  forward  elimination)  . Sim.ilarly,  back- 
ward recurrence  of  the  k-th  order  difference  equation 

' n 

starting  from  the  terminal  values  (which  are  all  zero)  produces  the  corre- 
sponding solution  of  the  boundary  value  problem;  this  is  the  back  substi- 
tution stage.  Data  is  given  in  the  tables  only  for  selected  values  of  r . 
All  of  the  values  of  r that  were  used  beyond  r = m are  included,  how- 
ever. These  values  illustrate  the  termination  procedure. 

It  should  be  noted  that  when  j = 4 the  method  as  implemented 
by  the  FORTRAN  program  is  equivalent  to,  but  not  identical  with, 
forward  recurrence  of  the  given  difference  equation.  See,  for 
example,  Table  8.  This  is  because  DGBFA,  the  UNPACK  matrix 
factorization  subroutine,  scales  the  linear  system  in  such  a way 


chat  the  principal  diagonal  of  the  lower  triangular  factor  is  a sequence 
of  ones.  Therefore,  the  result  of  the  forward  elimination  stage  differs 
from  the  result  of  the  back  substitution  stage  (unless  the  trailing  coeffic- 
ients of  the  original  linear  difference  operator  are  all  equal  to  one) . 

Since  there  is  no  question  of  truncation  error  when  j = 4 , the  sequences 

w (r)  , y (m)  and  ri  (m)  , r = m+l,m+2,...,  are  not  defined, 
m r ^ > i 

It  has  been  observed  that  Giver's  original  second-order  method  is 
likely  to  cause  overflow  in  most  computers.  This  same  likelihood  of  over- 
flow exists  for  the  extended  algorithm.  The  examples  to  be  presented  in 
this  chapter  were  computed  using  double  precision  on  a Univac  1108.  Double 

precision  was  used  not  for  higher  precision  but  for  wider  exponent  range. 

308 

The  double  precision  range  extends  to  10  whereas  the  largest  number 

200 

generated  in  the  examples  was  somewhat  greater  than  10  . Van  der  Cruyssen 

[25]  gives  a rescaled  version  of  the  second-order  algorithm,  based  on  LU 
decomposition,  that  greatly  reduces  the  incidence  of  overflow  or  underflow. 
An  equivalent  procedure  for  overcoming  this  difficulty,  based  directly  on 
the  use  of  ratios  of  consecutive  terms  in  each  solution,  is  described  by 
Aggarwal  and  Burgmeier  [1].  In  the  present  writer's  view,  however,  it  is 
preferable  (if  possible)  to  avoid  mathematical  modifications  of  algorithms 
simply  to  suit  arbitrary  limitations  of  existing  computers.  A software 
alternative  to  rescaling  is  proposed  in  [24],  in  whicn  a full  integer 
storage  location  is  allocated  to  the  exponent  of  each  floating-point 
number.  A FORTRAN  implementation  of  this  proposal  is  available;  see  [11]. 
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4 . 1 Examples  Involving  the  Homogeneous  JYIK  Equation 
From  Che  asymptotic  forms 


(4.1.1) 


J^(x) 


(2Trr) 


1/2 


\(x)  - 


and 


(4.1.2)  I^(x) 


(-)’^K^(x)  - (-) 


which  are  valid  for  fixed  x and  large  r,  we  obtain 

(4.1.3)  < I^(Xj)  < (-)'^K^(x,)  < Y^(x^) 
when  x^  < x^  and 

(4.1.4)  I^(xp  < J^(x^)  < Y^(x^)  < (-)''k^(x2) 

^jhen  x^  < , where  in  (4.1.3)  and  (4.1.4)  the  symbol  < is  used  to  indicate 

the  order  of  dominance  of  the  functions  as  r . Thus  when  x^  ^ x^ 

the  JYIK  operator  is  totally  separable.  As  with  any  totally  separable 
operator  of  order  4,  there  are  four  distinct  types  of  solution  and  the  dimen- 
sion of  the  subdominant  subspace  of  type  s is  s . This  means  that  j = s 
is  the  proper  number  of  initial  values  to  specify  for  a solution  of  type  s . 

When  all  four  Bessel  functions  have  the  same  argument,  from  (4.1.1) 
and  (4.1.2)  we  have 

(4.1.5)  J^(x)  ~ I^(x)  , (-)"K^(x)  --  (-)^^^  ^trY^(x)  , r ->  » . 

The  first  of  these  relations  suggests  that  we  define  the  further  solution 


Then  it  may  be  verified  that 


(4.1.7) 


A^(x) 


1 x~ 

2 r+1 


I^(x) 


J^(x) 


^ 00 


where  either  function  in  the  braces  on  the  right  may  be  chosen.  Further- 

more,  since  there  is  no  nonzero  linear  combination  of  (-)  K (x)  and 

r 

Y^(x)  that  reduces  to  a solution  of  lower  type,  we  conclude  that  there 
are  exactly  three  distinct  types  of  solution.  Symbolically,  we  have 


(4.1.8) 


r ^ 

I (x) 

^ r 

(-)  K (x) 

r 

< 

>•  < 

r 

> 

J (x) 

Y (x) 

L r J 

L ^ J 

where  in  each  pair  of  braces  either  function  may  be  chosen.  Thus  when 
x^  = x^  the  JYIK  operator  is  separable  but  not  totally  separable.  The 
dimensions  of  the  subdominant  subspaces  of  types  1,2  and  3 are  1,2  and 
4 respectively. 

Although  A^(x)  has  type  1,  it  is  separated  only  by  a relative 

factor  of  1/r  from  solutions  of  type  2.  Let  v (r)  denote  the  approxima- 

' n 

tion  to  A^(x)  that  is  obtained  from  the  algorithm  with  terminal  point  n. 
Then,  suppressing  the  arguments  of  A^(x)  and  all  the  Bessel  functions 
and  using  equation  (3.2.6),  we  have' 

y^(r)  ~ A^  - ?^(n)I^  - P2('^)Y^  ~ , n - » 

where 

Since  I^(x)  , Y^(x)  and  K^(x)  are  not  normalized  as  required  for  the 
validity  of  (3.2.6),  we  have  only  the  weaker  statement  given  here. 
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A 

n 

Y 

n 

/ 

I 

n 

Y 

n 

= 

"^n+1 

"^n+1 

Vi 

Vi 

^n+2 

/ 

n+2 

''n+2 

I 

n 

A 

n 

{-)\ 

/ 

I 

n 

Y 

n 

(-)\ 

"^n+1 

^nh2 

Vi 

^n+2 

^n+2 

/ 

^n+2 

and 


I 

n 

Y 

n 

A 

n 

/ 

I 

n 

Y 

n 

(-)\ 

p3(n)  = 

"^n+1 

/ 

T 

""n+l 

^n+2 

V2 

\+2 

/ 

^n+2 

V2 

In  consequence  of  (4.1.1),  (4.1.2),  and  (4.1.7)  we  see  that 
2 -1 

P^(n)  -X  ( 2n)  for  large  n , whereas  ^2^^^  ^3^^^  approach  zero 

much  more  rapidly.  Thus  we  expect  that  the  rate  of  convergence  of  the 
algorithm  will  be  approximately  as  1/n  . This  has  been  confirmed  in 
numerical  tests.  Since  the  separation  of  A^(x)  from  solutions  of  type 
two  is  so  weak,  it  is  usually  advantageous  to  regard  A^(x)  as  belonging 
to  the  same  type  as  d^(x)  and  I^(x)  ; this  means  that  we  specify  two 
initial  conditions  rather  than  one.  This  modification  improves  convergence; 
in  consequence  the  value  of  the  terminal  point  n is  reduced  very  consider- 
ably. This  gain  is  offset  to  a minor  extent  by  some  loss  in  accuracy  in 
the  values  of  A^(x)  for  higher  values  of  r ; the  extent  of  this  loss  is 
no  more  than  the  cancellation  that  takes  place  between  the  corresponding 

values  of  I (x)  and  J (x)  in  Eq.  (4.1,5). 
r r 


1 


Example  . 1 . 1 In  this  example  we  compute  the  solutions  J^(l)  , 1^(10)  , 
(-)'K^(IO)  and  Y^(l)  of  the  homogeneous  equation  defined  in  Figure  6, 
with  x^  = 1 , x^  = 10  , for  r = 0 , 1 , 2 , . . . , 100  . Since  each  of  these 
solutions  is  monotonic  (in  magnitude),  the  algorithm  will  be  numerically 
stable  provided  that  the  proper  number  j of  initial  values  is  used  in 
each  case.  This  number  is  j = 1 for  J^(l)  , j = 2 for  1^(10)  , 
j = 3 for  (-)^K^(IO)  and  j = 4 for  Y^(l)  ; compare  (4.1.3). 

Table  1 gives  the  numerical  coefficients  of  the  JYIK  operator  with 
X.,  = 1 and  = 10  for  r = 0 , 1 , 2 , . . . , 109 . Tables  2,  3 and  4 give  the 
numerical  coefficients  of  the  (0 ,j ) -factorizations  of  this  JYIK  operator 
for  j =1,2,  and  3.  (We  do  not  present  the  ( 0, 4) -factorization . ) Since 

X X x+i 

the  (i, j)-factorization  = A^B  produced  by  the  test  program  has 

the  principal  diagonal  of  equal  to  (1,1,1,...)  , we  do  not  include 

values  of  a^ (r)  . 

Tables  5 through  8 give  the  results  of  the  computation  of  the  four  ' 

desired  solutions  J^(l)  , 1^(10)  , (-)^K^(IO)  and  Y^(l)  for  selected 
values  of  r ; see  §4.0  for  a full  description  of  the  tables.  Initial 
%»alues  were  taken  from  Tables  9.4  and  9.11  in  [16];  they  are  given  in 
parentheses.  These  tables  also  supply  values  of  the  four  Bessel  functions 
for  certain  other  values  of  r in  the  range  0 5 r 5 100  . These  entries 
all  agree  with  our  computed  results,  except  for  the  tenth  significant 
figure  in  some  instances.  This  slight  discrepancy  is  explained  by  the 
fact  that  the  starting  values  that  were  used  inevitably  have  an  error  of 
several  units  in  the  eleventh  significant  figure.  Incidentally,  because  in 
this  (and  other  examples)  we  have  not  attempted  rigorous  control  of  I 

I 

) 

( 


individual  rounding  errors,  a guaranteed  error  bound  cannot  be  given. 


Nevertheless,  the  stability  of  the  algorithm  for  all  solutions  in  this 


example  is  firmly  demonstrated. 
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Table  4.  Numerical  coefficients  of  tlie  (0, 3)-factorization  of  the  JYIK  operator  with  x = 1 and  x - 10. 
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Table  6.  Computation  of  I (10)  from  the  JYIK  equation  with  x -1,  x -10,  1-0,  m-100 
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Table  7.  Computation  of  (-)  K (10)  from  the  JYIK  equation  with  x =1,  x =10,  1=0,  m=100. 
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Examp 1 e 4.1.2  In  this  example  we  compute  the  solutions  1^(1)  , J (10)  . 
Y^(IO)  and  (-)  K^(l)  of  the  homogeneous  equation  defined  in  Figure  6, 
with  x^  = 10  , x^  = 1 , for  r = 0 , 1 , 2 , . . . , 100  . Both  1^.(1) 

K^(l)  are  monotonic  over  this  range  of  r . The  solutions  J (10) 
and  Y^(IO)  are  monotonic  only  for  10  < r 5 100  ; for  lower 
values  they  oscillate  with  approximately  the  same  amplitude.  However, 
since  there  is  no  "crossing  over"  of  growth  behavior  in  the  four  solutions 
in  the  nonasymptotic  region  0 < r < 10  , we  expect  the  algorithm  to  be 
stable  with  the  following  numbers  of  initial  conditions:  j = 1 for  1^(1), 
j = 2 for  J^(IO)  , j = 3 for  Y^(IO)  and  j = 4 for  (-)’^K^(l)  ; 
compare  (4.1.4).  Tables  9 through  12  give  the  results  of  the  computation 
of  each  of  these  solutions  for  selected  values  of  r ; see  §4.0  for 
a full  description  of  these  tables.  Again,  the  stability  of  the  algorithm 
for  each  solution  is  confirmed  by  comparison  with  Tables  9.4  and  9.11  of 


[16]. 
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Table  9.  Computation  of  I (1)  from  the  JYIK  equation  with  x =10,  x =1,  i=0,  m=100, 
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Table  10.  Computation  of  J (10)  from  the  JYIK  equation  with  x =10,  x„=l,  1=0,  m=100 
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Example  4.1.3  In  this  example  we  compute  the  solutions^ A^(l)  , J^(l)  , 
1^(1)  , Y^(l)  and  (-)  K^(l)  of  the  homogeneous  equation  defined  in 
Figure  6,  with  x^  = 1 , x^  = 1 , for  r = 0 , 1 , 2, . . . , 100  . Since  each 
solution  is  monotonic  (in  magnitude)  except  possibly  at  the  very  beginning 
of  the  range  of  r , the  algorithm  is  expected  to  be  stable  provided  that 
the  proper  number  of  initial  values  is  taken  in  each  case.  This  number 
is  j = 1 for  A^(l)  , j = 2 for  J^(l)  and  ^ ~ ^ f°^ 

X 

Y^(l)  and  (-)  K^(l)  • However,  in  the  case  of  A^(l) , for  reasons  given 
earlier  in  this  section,  we  take  j = 2 instead  of  j = 1 , thereby 
exchanging  mild  instability  for  a much  improved  rate  of  convergence. 

Because  of  this  mild  instability,  each  individual  rounding  error, 
introduced  at  say  r = s,  is  amplified  in  direct  proportion  to  the  number 
of  steps  taken  beyond  the  point  s , since  in  (4.1,7)  A^/ = 0(l/r)  . 

Thus,  after  100  steps,  approximately  two  decimal  places  of  precision  will 
be  lost.  On  the  other  hand,  if  we  were  to  take  the  number  of  initial 
conditions  theoretically  required  (in  this  case  it  is  J~=l,-  since 
is  dominated  by  all  of  the  other  homogeneous  solutions),  then  the  conver- 

g 

gence  would  be  extremely  slow;  A terminal  point  of  the  order  of  10 
would  be  required  to  obtain  8 significant  figures  of  accuracy  in  A^^^d)  . 

Tables  13  through  17  give  the  results  of  these  computations  for  selected 
values  of  r ; see  §4.0  for  a full  description  of  these  tables.  Again, 
these  results  agree  with  the  entries  in  Table  9.4  and  9.11  of  [16].  In 
the  case  of  A^(l)  , the  method  we  have  used  is  equivalent  to  computing 
J (1)  and  I (1)  as  type  2 solutions  and  then  forming  their  difference. 


The  function  A (1)  is  defined  by  A (1)  = J (1)  - I (1)  ; compare 
(4.1.6).  ^ r r r 
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Table  14.  Computation  of  J (1)  from  the  JYIK  equation  with  x =x  =1,  i-0,  m=100. 
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4 . 2 Examples  Involving  the  Inhomogeneous  JYIK  Equation 

For  fixed  x and  large  r the  Anger-Weber  function  has  the 
asymptotic  forms 


(4.2.1) 


EzrCx) 


2x 


(4r  -1)TT 


compare  [21].  Comparison  of  (4.2.1)  with  (4.1.1)  and  (4.1.2)  shows  that 


'lyxz)' 

1 

(4.2.2)  . 

ycxy 

J 

► < E^(x)  < .< 

► 

vyxi) 

for  any  given  set  of  values  of  x^  , x^  , x . Thus  E^(x)  is  a solution 
of  type  2 of  the  first  inhomogeneous  equation  defined  in  Figure  6 and  the 
correct  number  of  initial  values  of  E^(x)  to  be  used  in  every  case  is 
j = 2 . 

The  asymptotic  form  of  the  Struve  function,  also  taken  from  [21],  is 


(4.2.3) 


y 2 ' ' 


rrr 


for  fixed  x and  large  r . Comparison  of  (4.2.3)  with  (4.1.1)  and  (4.1.2) 


shows  that 
(4.2.4) 


ry(x)'i 


H (x) 


Vm 


> , r 


L 'J^(x)  J 


Thus  H^(x)  is  a solution  of  type  1 of  the  second  inhomogeneous  equation 
defined  in  Figure  6 when  = x^  ; compare  (4.1.8).  However,  since 

the  separation  of  H^(x)  from  solutions  of  type  2 is  even  weaker  than 


chac  of  the  solution  A^(x)  defined  by  (4.1.6)  it  is  again  advantageous 
in  most  applications  to  specify  the  number  of  initial  values  appropriate 
to  solutions  of  type  2,  i.e.,  j = 2 instead  of  j = 1 ; compare 
Example  4.1.3. 

Example  4.2.1  For  the  inhomogeneous  fourth-order  difference  equations 
defined  in  Figure  6,  we  compute  E^(l)  with  x^  = x^  = 1 for  r = 0,1,2, 
...,100  ; and  H^(O.l)  with  ~ ^2  ~ ^ ~ 0,1 , 2,  . . . , 50 . ' 

The  appropriate  number  of  initial  values  in  these  cases  is  j = 2 for 
E (1)  and  j = 1 for  (0.1)  . However,  for  reasons  given  above  we 

actually  use  j = 2 for  H^(O.l).  Initial  values  for  E^(l)  and  H^(O.l) 

were  computed  from  series  expansions  given  in  [21]  and  checked  against  value 
found  in  [21]  and  [23].  Tables  18  and  19  give  the  results  of  these  com- 
putations for  selected  values  of  r ; see  §4.0  for  a full  description  of 
these  tables.  They  agree  with  the  10-figure  values  of  E^(l)  and  H^(O.l) 
r = 2,3,4, 5 , given  in  [23];  the  7 through  9-figure  values  of  E^(l)  , 
r = 2, 3,..., 10  , given  in  [21];  and  the  9-figure  values  of  H^(O.l)  , 
r = 1,2,..., 13  , given  in  [21].  In  addition,  they  agree  with  series  calcu- 
lation of  the  functions  to  10  significant  figures  in  the  case  of  E^(l)  , 
r = 2, 3,..., 100  , and  to  9 significant  figures  in  the  case  of  H^(O.l)  , 
r = 2, 3,..., 50  . The  loss  of  precision  in  the  Struve  function  computa- 
tion IS  due  to  mild  instability  which  is  acceptable  in  lieu  of  impossibly 
slow  convergence. 


We  restrict  the  sequence  of  H (0.1)  to  r < 50  in  order  to  avoid  under- 
flow in  the  computer.  As  suggested  at  the  end  of  §4.0,  this  difficulty 
could  be  overcome  by  using  the  software  described  in  [11]. 
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Table  18.  Computation  of  £ (1)  from  the  JYIK  equation  with  x =x  =1,  i=0,  m=100 
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4 . 3 Examples  Involving  the  HomogeneotiS  JYPQ  Equation 

For  large  values  of  r and  fixed  values  of  p and  x , the  associated 
Legendre  functions  have  the  asymptotic  forms 

(4.3.1)  P^(x)  = ^ ^^^cos(r0  + ^ - ^ + + O(^) 


and 


(4.3.2) 


Q!^(x) 


r(r^+l) 

r(r+3/2) 


(■ 


TT_ 

2 sin  G 


1/2 


■) 


cos (r 0 + —{ 


+ ^ + iin) 


+ 0(i) 


where  x = cos  0 with  0 < 0 < rr  ; see  [16,  Chapter  8]  . Consequently, 
for  values  of  x less  than  1 in  magnitude,  neither  associated  Legendre 
function  dominates  the  other.  Comparing  (4.3.1)  and  (4.3.2)  with  (4.1.1) 
we  obtain 


"P^(X2)" 


(4.3.3) 


Jr(Xi)  < 


Q^(k2) 


► < » r ->  « 


for  all  admissible  values  of  [j,  , x^  , X2  • Thus  the  JYPQ  equation  has 
three  distinct  types  of  solution,  and  the  dimensions  of  the  subdominant 
subspaces  of  types  1,  2 and  3 are  1 , 3 and  4 respectively;  this 
contrasts  with  (4.1.3),  (4.1.4),  and  (4.1.8).  ■ 

Example  4.3.1  In  this  example  we  compute  the  solutions  3^(1)  ? P^(0.5)  , 
Q^(0.5)  and  Y^(l)  of  the  homogeneous  equation  defined  in  Figure  7, 
with  u=0  , x^  = l , X2  = 0.5  , for  r = 0 , 1 , 2 , . . . , 100  . Since  these 
solutions  exhibit  their  characteristic  asymptotic  separation,  described  by 
(4.3.3),  from  the  very  beginning  of  this  range  of  r , the  algorithm 


"'"in  accordance  x^ith  custom,  when  the  order  of  the  associated  Legendre 
functions  is  zero  we  drop  the  superscript  in  the  notation. 


will  be  stable  provided  the  proper  number  of  initial  values  is  used  in 

each  case.  This  number  is  j = 1 for  d^(l)  , j = 3 for  P (0.5)  and 

Q^(0.5)  , and  j = 4 for  Y^(l)  . 

Table  20  gives  the  numerical  coefficients  of  the  JYPQ  operator  with 

u = 0 , = 1 and  = 0.5  for  r = 0 , 1 , 2 , . . . , 104 . Tables  21  and  22 

give  the  numerical  coefficients  of  the  (0 , 1) -factorization  and  the  (0,3)- 
factorization . Since  the  (i,  j )-factorization  produced  by 

the  test  program  has  the  principal  diagonal  of  equal  to  (1,1,1,...), 

we  do  not  include  values  of  a^ (r)  . 

Tables  23  through  26  give  the  results  of  the  computation  of  the  four 
desired  solutions  for  selected  values  of  r ; see  §4.0  for  a full  descrip- 
tion of  these  tables.  Initial  values  were  obtained  from  Table  9.4  of  [16] 
in  the  case  of  the  Bessel  functions  and  by  hand  calculation  of  the  explicit 
elementary  forms  given  in  [16],  Section  8.4  in  the  case  of  the  Legendre 
functions.  The  computed  Bessel  function  values  agree  satisfactorily  with 
the  values  given  in  [16] . The  Legendre  function  values  were  compared  with 
data  available  for  r < 10  in  [16]  and  (in  the  case  of  the  P^(0.5))  [17]. 
Agreonent  is  satisfactory  to  the  precision  available  in  these  references; 
this  ranges  from  6 to  8 significant  figures.  In  addition,  all  of  the 
computed  Legendre  function  values  agree  satisfactorily  with  double-preci- 
sion values  computed  by  an  unpublished  program  of  Dr.  J.  M.  Smith  of  the 
National  Bureau  of  Standards.  This  is  sufficient  to  demonstrate  stability. 
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Table  21.  Numerical  coefficients  of  the  (0,  l)-f  actorizatlon  of  the  JYPQ  operator  with  p=(),  x =1  and 
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Table  22,  Numerical  coefficients  of  the  (0, 3)-factorization  of  the  JYPQ  operator  with  p=0,  x -1  and 
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Table  23.  Computation  of  J (1)  from  the  JYPQ  equation  with  p=0,  x =1,  x =0.5,  i=0,  m=100 
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able  24.  Computation  of  P (0.5)  from  the  JYPQ  equation  witli  p=0,  x^=l,  x =0.5,  1=0,  m=100 
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Table  25.  Computation  of  Q (0.5)  from  the  JYPQ  equation  with  p=0,  x =1,  x„=0.5,  1=0,  m=100 
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Example  4.3.2 


In  this  example  we  compute  the  solutions  P^(0.5)  and 
Q^(0.5)  of  the  homogeneous  equation  defined  in  Figure  7,  with  la  = 1 , 
x^  = 1 , = 0.5  , for  r = 0 , 1 , 2 , . . . , 100  . The  algorithm  will  be  stable 

provided  that  we  take  three  initial  values  for  each  of  these  functions; 
compare  (4.3.3).  Tables  27  and  28  give  the  results  of  these  computations 
for  selected  values  of  r ; see  §4.0  for  a full  description  of  these 
tables.  Initial  values  were  obtained  by  hand  calculation  of  the  explicit 
elementary  forms  given  in  Tables  8.2  and  8.4  of  [16].  The  computed  function 
values  agreed  satisfactorily  for  r 5 10  with  the  limited  data  available 
in  [16]  and  [17] . In  addition,  satisfactory  agreement  was  obtained  with 
double-precision  values  produced  by  an  unpublished  program  of  Dr.  J.  M. 
Smith;  compare  Example  4.3.1.  Consequently,  the  stability  of  the  method 


is  confirmed  for  this  example. 
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4 . 4 Examples  Involving  the  Inhomogeneous  JYPQ  Equation 
Using  Stirling's  formula,  we  find 


(4.4.1) 


r(r+p+l) 

r(r+3/2) 


U-1/2 


for  fixed  values  of  u = 0,1,2,...  . Comparison  of  (4.2.1)  with  (4.3.1), 
(4.3.2)  and  (4.4.1)  shows  that 


(4.4.2) 


Er(xi)  < ^ 


Q^(x2) 


for  fixed  [X  , x^  and  x^  such  that  x^  > 0 and  ^9  < ^ • Accordingly, 

the  Anger-Weber  function  is  a solution  of  type  1 and  asymptotically  the 

correct  number  of  initial  values  to  be  used  in  every  case  is  j = 1 ; 

compared  (4.3.3).  However,  the  separation  of  E^(x^)  from  the  solutions 

-1/2 

of  type  2 is  weak  when  u = 0 (being  of  relative  order  r at  the  odd 

-3/2 

terms  of  the  sequence  and  r at  the  even  terms  of  the  sequence) , The 

-3/2 

separation  is  also  weak  when  u-  = 1 (being  of  relative  order  r at 

-5/2 

the  odd  terms  of  the  sequence  and  r at  the  even  terms  of  the  sequence) 

Accordingly,  we  treat  E^(l)  as  a solution  of  type  2 in  both  cases,  even 
though  when  u,  = 1 greater  numerical  instability  is  incurred  than  in  any 
of  the  other  examples.  This  is  especially  true  at  the  even  terms  of  the 

sequence,  which  are  separated  from  the  Legendre  functions  by  the  relative 

a P -5/2 
order  of  r 

From  (4.2.4)  and  (4.3.3)  it  is  seen  that  the  type  of  the  Struve 
function  H (x^)  with  respect  to  the  JYPQ  operator  is  0 for  all 
admissible  values  of  x^  and  x^  . But  again  because  of  the  weakness 
of  the  separation,  we  regard  H^(x^)  as  belonging  to  type  1. 


Example  4.4,1  For  the  inhomogeneous  fourth-order  difference  equations 
defined  in  Figure  7,  we  compute  E^(l)  with  u = 1 , = 1 and  x^  = 0.5  , 

for  r = 0,1,2,  ...  ,100;  and  H^(O.l)  with  u.  = 1 , x^  = 0.1  and  x^  = 0.5 
for  r = 0,1,2, ... ,50.  (We  restrict  the  sequence  of  H^(O.l)  to  r S 50 

for  the  same  reasons  that  were  described  in  Example  4.2.1.)  For  reasons 
given  above,  we  use  3 initial  values  (rather  than  1)  for  E^(l)  and  1 

initial  value  (rather  than  0)  for  H (0.1)  . Initial  values  for  E (1) 

r r 

and  H^(O.l)  were  computed  from  series  expansions  given  in  [21]  and 

checked  against  values  found  in  [21]  and  [23]. 

Tables  29  and  30  give  extracts  from  the  results  of  these  computations 

for  selected  values  of  r ; see  §4.0  for  a full  description  of  these 

tables.  Comparison  of  the  function  values  produced  by  the  algorithm  with 

values  given  in  [21]  and  [23]  show  good  agreement  over  the  available 

ranges,  which,  however,  do  not  extend  beyond  r = 13  . Comparison  of  the 

values  of  E^(l)  in  Table  2.9  with  values  computed  to  10  significant 

figures  from  the  series  expansion  shows  a steady  deterioration  of  precision 

at  the  even  values  of  r , culminating  in  a relative  error  of  0.5(10)  ^ 

at  r = 100  . At  the  odd  values  of  r the  precision  also  deteriorates 

steadily,  but  at  a slower  rate,  culminating  in  a relative  error  of 

0.8(10)  ^ at  r = 99.  These  rates  of  loss  of  precision  are  consistent 

-5/2  -3/2 

with  the  separation  ratios  of  r and  r which  govern  the  growth 


of  the  error  in  this  example. 


-10 


That  is,  an  individual  relative  rounding  error  of,  say,  10  in  one 

of  the  stored  initial  values  of  E (1)  will  be  amplified  by  a factor  of 

S/75  3/23 

approximately  100  = 10  at  E _(1)  and  by  100  = 10  at  E (1). 

i uu  y y 

(As  noted  in  Eq.  (4.2.1),  the  function  is  in  behavior  different  at  even  and 
odd  values  of  r .)  These  amplifications  after  100  steps  are  greater  than 
the  factor  of  100  that  was  predicted  and  observed  in  the  computation  of 
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A (1)  in  Example  4.1.3  and  Table  13,  where  A (1)  was  computed  as  if 
r r 

it  were  a function  of  the  next  higher  type  than  it  actually  is. 

Finally,  comparison  of  the  values  of  H^(O.l)  in  Table  30  with 

values  computed  to  10  significant  figures  from  the  series  expansion 

shows  agreement  to  9 or  more  significant  figures  at  all  values  of  r . 

The  loss  of  precision  in  the  Struve  function  computation  is  slight,  due 

-1/2 

to  the  separation  ratio  of  only  r , but  it  offsets  impossibly  slow 

convergence  which  would  be  evidenced  if  the  theoretically  correct  number 
of  initial  values  were  to  be  taken. 


o 

— 

— 

M 

'v 

to 

A 

V 

5^' 

£0 

O 

'N 

N.* 

C i 
O 

</ 

/K1 

V 

jM 

CJ 

'J 

V V 

>V 

<►  y 

9 

'O 

o 

Q 

9 

C^‘) 

o 

9 

o 

9 

Q 

9 

iJ> 

9 

9 

1 * 

•N 

/•  ■^ 

o 

o 

AT 

/\ 

N 

— 

-9 

■j'. 

,-s 

•Nj 

9. 

J.. 

> 

— 

V 

d) 

« 

-V 

'V 

,\» 

o 

Aj 

1 .) 

‘s'i 

■*--m 

^si 

o 

:-J 

r'; 

O 

Txs 

r ^ 

C] 

> 

VJ 

-J 

J 

^7' 

■/^ 

;-J 

C»5 

;*> 

7> 

O' 

• > 

r; 

7^* 

-■N 

ii* 

O 

.^i 

‘s/ 

'o 

- A 

«a4 

s . 

•r-4 

C-) 

n.-* 

V 

7o 

" ] 

: ^ c'7 

■» 

V 

— 

■J 

-:T- 

-OJ 

Uv 

” 

V.? 

C'l 

o 

^‘3* 

W 

■'J 

o 

ffi 

’^7’’ 

« 

.>.A 

'_X> 

o 

N 

■I' t 

'■' ';  w 

L'^ 

</ 

Cv 

-7' 

O': 

C' 

C> 

ft) 

v<? 

\Q 

'T.i 

>;> 

Q'~) 

\5 

•»K 

2'J 

-■  N 
N.  V 

•— < 

“^v 

• 

Cw 

■V 

o 

Vv 

N.>* 

C-^ 

^0 

vO 

,'] 

V 

.0 

Q 

o 

/»  A 

CO 

3 

<^•7 

>»*3> 

•7  ;“*? 

m 

V 

N 

c: 

« 

iO 

N 

(■'.j 

a 

t'i 

N 

S' 

V* 

CO 

CO 

>r- 

ro=» 

o o 

II 

AT 

V 

-I.' 

I'’ 

Vi 

* 

;s_. 

> 

>■» 

''iP 

7-» 

o 

O 

•.P 

■'i-5 

o 

\w 

•7. 

■•>  C.5 

•r-) 

<? 

S' 

f?\ 

G C3 

O 

o 

o 

— 

■— . ^ 

/<— V 

G 

1 

0 

1 

Y 

0 

1 

© O 
1 1 

o 

G 

<? 

'J' 

• A 

O (S' 

IS 

03 

LS 

o 

G N 

C? 

uo 

G 

tD 

© 

G 

S' 

—1 

CO  CO 

iM 

N. 

© 

ClJ 

O !S. 

03 

o 

•c? 

T'j  ?0 

S' 

© 

o 

:\.  N. 

LS 

o 

L'tl 

- LS 

o 

to 

N 

03 

« 

N 

-< 

« CO 

T? 

■S' 

S' 

s 

S' 

c 

o 

G 

© 

o 

© 

o 

0 

1 

0 

1 

9 

O 

1 

9 

9 

0 

1 

/«— V 

i 

lit) 

1 

L=5 

I 

9 

1 

o 

( 

I 

CO 

1 

CO 

o 

03 

i'* 

t'* 

C3 

C3 

o 

»-* 

S' 

© 

e. 

© 

S' 

N 

N 

G 

(03 

C3 

C3 

C-3 

'w' 

o 

CO 

O 

CO 

C'O 

CO 

C'O 

'C' 

C3 

S' 

S' 

S' 

S' 

>T 

Ifl 

N 

'm 

N 

3'* 

is 

N 

''P 

© 

'v 

G 

© 

'3 

03 

Cv 

CO 

(© 

CO 

CO 

Cl 

C3 

'3 

'S 

o 

'3 

v3 

(Cl 

'<5 

o 

o 

O 

o 

c 

o 

O 

o 

o 

o 

© 

9 

vO 

'0 

o 

CJ 

6 

LO 

N 

N 

'3 

LO 

^7"* 

■N 

U 

LO 

VO 

LS 

LO 

N 

LO 

© 

C'] 

CO 

V 

SI 

s> 

S' 

S' 

V 

LO 

C'J 

o 

N 

LO 

Ni 

CO 

© 

wmM 

o 

N 

t' 

—1 

O 

VO 

\~4 

S' 

(M 

V 

Cl 

S' 

c c 

© 

Cl 

!■? 

O 

LO 

.- 

Vj 

N, 

<•«% 

A"* 

C3 

W-J 

A.« 

o 

— O 

S-H  Cv 

o 

■-  o 

© 

— 

C' 

© 

o o 

y\ 

o 

o 

o o 

77 

© 

Q 

b 

o o 

AS 

O 

o o 

o 

o 

’-'  o 

— CO 

T— * 

CJ 

’-'  o 

»— 

o 

o 

CJ 

o 

o © 

o 

o 

.rs 

✓ 'S 

G 

c 

c 

CO 

o 

G O 

o 

o o 

o o o 

O CO 

^ 

o 

COG 

o 

CJ 

G 

o 

o 

-)- 

.'• 

-J- 

-'- 

J. 

“S* 

+ 

•r 

-,- 

V 

-,- 

•• 

*r  *r 

T *•* 

-I- 

"r  "* 

-,- 

**” 

•,“ 

C Cl 

n 

r >4 

,ST 

c>  Cl 

o 

ct 

' ^ 

•■O 

CD 

CO 

*''7^ 

O C3 

S'  C -1 

CO 

'■.0  CO 

S' 

CD 

6 

t-  --o 

Cl 

'.“t 

o 

r »• 

Cl 

© 

C3 

t-- 

c> 

o 

C7 

o 

C'l  N 3- 

Cl  Cl  O O 

o 

O N 

C' 

’>■? 

-J 

CO 

C-'' 

S' 

/■  7 ' • 

•A 

'J' 

w 

\0 

\. 

01 

LO 

A*, 

Vv 

A'. 

© ©■ 

LO 

C'j 

o 

LO 

Cl 

tj' 

N 

C "r 

'"; 

; * 

y s 

As 

CD 

C'» 

C'O 

A> 

o 

c* 

. J 

-^1 

vC  O'* 

O CO 

o 

O -1 

^/■' 

N 

CO 

■c '.: 

— 

y •< 

r*4 

^ S 

r '> 

«« 

'■  1 

O t-' 

— • 

• ' 

C'O 

o 

A** 

» A 

S' 

S'  C 3 vO  L-3 

L3 

\* 

CD  Cl 

a 

S' 

J*  'l 

■p- 

u 

L'.(  0 

As 

r 

'7* 

C** 

r • 

7 

l; 

Cj 

f/- 

Lt 

t ■' 

“,  J 

J-*; 

p.  A-  ] 

VO  •' 

J,' 

C' 

C-j 

L.' 

O 

'i  I 

•M 

c 

'T*. 

ijj 

fn. 

•\0 

|XS 

t J 

•».1 

■> 

S' 

C5 

C'  V 

<01 

to 

o 

o -■' 

CO 

N 

v-'j 

^0 

L'; 

r 1 '•'' 

'A 

2 J 

7^ 

••A. 

CO 

cc 

CO 

CD 

w' 

k-J 

7o 

-'  C i 

in 

CJ 

CO 

LO  Cl 

CJ 

\’J 

C.o 

C’O 

VI) 

O 

O CO 

i •» 

.*• 

Vl  C'1 

V 

Y'xi 

r ^ 

•’J-' 

o 

i'-rt 

LO 

N 

C'O  Cj 

0 C,. 

wm 

c-l 

LO  S' 

o 

\\J 

LO 

CD 

t'i 

CO 

*“■ 

mm 

C-3 

fM 

1.0 

LO 

S' 

•“* 

Cl 

r 5 r«M 

CO  O'l 

<01 

S'  Cl 

LO 

Cl 

Cl 

'J 

CO 

1 

C — C!  C'O  S'  LO  3 N CD  S'  G © C © © O O CJ 

C - Cl  G S'  LO  -’J  N CD  © O 

— S’"  t:") 

— Cl  CO  S'  LO  O N LO 

S'  C'  S'  S'  S'  S'  S'  © S'  S'  o 

o c o o o 

c 

Table  29.  Computation  of  E (1)  from  the  JYPQ  equation  with  p=l,  x =1,  x =0.5,  1-0,  m-100 
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CHAPTER  5.  SUMMARY  AND  CONCLUSIONS 


This  thesis  is  concerned  with  the  numerical  solution  of  a broad  class 
of  linear  difference  equations  of  arbitrary  order,  and  of  either  homogeneous 
or  inhomogeneous  form.  An  algorithm  is  presented  for  the  stable  computa- 
tion of  any  solution  of  such  an  equation.  The  results  apply  directly  to 
scalar  linear  difference  equations  but  could  be  adapted  to  similar  problems 
posed  in  vector  form. 

Historically,  the  algorithm  is  an  outgrowth  of  algorithms  of  Miller 
and  Olver.  Tne  underlying  idea  of  Miller's  algorithm  is  the  following: 

If  a solution  y(r)  has  growth  behavior,  as  r ->■  ®,  that  is  less  than  or 
equal  to  the  growth  behavior  of  all  independent  solutions,  then  in  the 
reverse  direction  its  growth  behavior  is  greater  than  or  equal  to  the 
growth  behavior  of  all  independent  solutions.  Consequently,  because  of 
the  superposition  principle,  errors  introduced  in  computing  y(r)  by 
backward  recurrence  cannot  grow  faster  than  y(r)  itself,  and  in  this  sense 
the  backward  recurrence  process  is  stable.  However,  this  is  true  only  in 
a region  in  which  the  solutions  actually  maintain  the  growth  behavior 
indicated  by  their  asymptotic  forms  as  r For  this  reason  we  call 

this  the  asymptotic  region;  evidently  it  depends  on  the  asymptotic  struc- 
ture of  the  given  difference  equation  for  large  r. 

Miller's  algorithm  starts  the  backward  recurrence  process  at  a 
value  of  r that  is  larger  than  any  r for  which  a valid  approximation 
of  y(r)  is  wanted.  Arbitrary  starting  values  are  assumed.  The  computed 


145 


solution  is  essentially  a scalar  multiple  of  y(r)  for  early  values  of  r 
when  (i)  the  difference  equation  is  homogeneous;  (ii)  y(r)  is  subdominant 
compared  to  every  independent  solution  of  the  equation.  Under  these  con- 
ditions the  algorithm  is  stable  and  y(r)  is  obtained  from  the  computed 
solution  for  early  values  of  r by  a normalization  procedure,  such  as 
matching  the  value  of  the  computed  solution  at  r = 0 with  y(0)  in  order 
to  determine  the  correct  scale  factor.  Modified  normalization  procedures 
are  required  when  the  asymptotic  behaviors  of  more  than  one  fundamental 
solution  are  similar. 

Miller's  algorithm  involves  no  forward  recurrence,  and  is  stable 
only  for  solutions  whose  growth  behavior  is  less  than  or  equal  to  the 
growth  behavior  of  all  independent  solutions.  Giver's  algorithm  was 
applied  originally  to  second-order  inhomogeneous  difference  equations, 
in  cases  where  a particular  solution  y(r)  has  growth  behavior  as  r 
that  lies  between  that  of  a pair  of  linearly/  independent  solutions  of 
the  homogeneous  equation.  To  construct  a stable  algor ithm, Giver  posed 
a boundary  value  problem  having  one  condition  at  each  end.  The  initial 
condition  is  the  value  of  y(r)  at  a suitable  point  r = i.  The  terminal 
condition  is  zero,  it  being  anticipated  that  isolated  values  of  y(r) 
for  large  r will  not  be  easily  obtained.  The  resulting  boundary  value 
problem  is  a tri-diagonal  linear  system  of  algebraic  equations.  Giver 
specified  a method  for  solving  this  linear  system  that  involves  a 
forward  recurrence  stage  and  a bacfavard  recurrence  stage.  The  method 
is  equivalent  to  Gaussian  elimination  without  pivoting,  scaled  in  a 


particular  way.  Analysis  of  the  method  shows  that  (i)  the  forward  and 
backward  recurrence  stages  are  stable,  at  least  in  ranges  beginning  at 
a sufficiently  large  value  of  r;  (ii)  the  second-order  linear  difference 
operator  that  corresponds  to  the  difference  equation  is,  in  effect, 
factored  into  the  product  of  two  first-order  operators.  Furthermore 
Olver  derived  an  infinite  series  expansion  of  the  truncation  error,  valid 
at  any  point  r,  that  is  incurred  by  choosing  a particular  terminal 
point  for  the  boundary  value  problem.  This  enabled  the  optimal  terminal 
point  to  be  calculated  automatically,  in  contrast  to  Miller's  algorithm 
in  which  the  terminal  point  is  guessed  and  subsequently  tested. 

Extensions  of  Giver's  algorithm  to  difference  equations  of  higher 
order  were  proposed  by  Zahar,  Oliver  and  Cash.  Criteria  for  selecting 
the  correct  nimnber  of  initial  conditions  (and  therefore  the  correct 
number  of  terminal  conditions,  since  the  total  number  of  conditions 
must  equal  the  order  of  the  difference  equation)  were  given  independently 
by  Zahar  and  Oliver.  Zahar  proved  convergence  of  the  algorithm  under 
assumed  conditions  on  the  solutions  of  the  adjoint  equation.  Oliver 
investigated  stability  by  analyzing  the  factorization  of  the  finite 
boundary  value  problem.  Cash  provided  an  extension  of  the  series 
expansion  of  the  truncation  error.  The  work  of  Zahar,  Oliver  and  Cash, 
taken  together,  contains  the  essential  elements  of  a practicable 
generalization  of  Giver's  algorithm. 

In  this  thesis  Giver's  algorithm  has  been  extended  and  analyzed 
afresh.  The  view  was  taken  that  a linear  difference  operator  is  an 
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infinite  upper-triangular  band  matrix.  For  a specified  initial  point  i 
and  number  of  initial  conditions  j,  an  (i, j )-factorization  of  a linear 
difference  operator  D is  said  to  exist  if  Gaussian  elimination  without 
pivoting  about  the  j-th  super-diagonal  of  D,  starting  at  the  point  i 
on  this  diagonal,  can  be  carried  out  indefinitely.  This  process  produces 
two  lower-order  difference  operators  (infinite  band  matrices) . Their 
product  is,  of  course,  the  original  operator,  and  the  sum  of  their 
orders  is  the  order  of  the  original  operator.  Furthermore,  the  solution 
of  any  finite  boundary-value  problem  having  initial  point  i,  number  of 
initial  points  j and  terminal  point  n may  be  found  by  forward  recursion 
with  one  of  these  lower-order  operators  followed  backward  recursion  with 
the  other  (Theorem  1.2.2).  Procedurally,  a finite  boundary  value  problem 
is  solved  by  the  algebraic  method  of  forward  elimination  without  pivoting 
followed  by  back  substitution;  this  is  accompanied  by  an  analysis  of 
the  stability  of  this  process  in  terms  of  linear  difference  operators. 

Determination  of  the  correct  number  of  initial  conditions  for 
stability  requires  a knowledge  of  the  relative  rates  of  asymptotic 
growth  of  solutions  of  the  homogeneous  equation  as  r ^ . Since  the 

number  of  linearly  independent  homogeneous  solutions  is  equal  to  the 
order  il  of  the  linear  difference  operator,  there  cannot  be  more  than 
I distinct  growth  rates  as  r Cases  in  which  there  are  exactly 

I distinct  growth  rates  are  the  easiest  to  treat.  In  these  cases  a 
basis  exists  that  is  linearly  ordered  by  growth  rate,  and  we  call  the 
operator  totally  separable.  However,  there  exist  operators  of  practical 


interest, including  constant-coefficient  operators  for  which  one  or  more 
of  the  distinct  characteristic  roots  are  equal  in  magnitude,  that  are 
not  totally  separable.  Accordingly,  a more  general  class  of  operators, 
called  separable  operators,  is  introduced  (Definition  2.2.1).  The  solu- 
tions of  a homogeneous  equation  that  has  a separable,  but  not  totally 
separable,  operator  exhibit  less  than  Z distinct  rates  of  growth  as  r ^ “ 
Separable  operators  are  characterized  by  the  existence  of  certain  bases, 
called  optimally  ranked.  An  optimally  ranked  basis  includes  a ^ I 
disjoint  subsets  such  that  (i)  neither  of  two  distinct  solutions  in 
the  same  subset  dominates  the  other  as  r (ii)  any  two  solutions 

taken  from  distinct  subsets  are  such  that  one  dominates  the  other; 

(iii)  no  nonzero  linear  combination  of  solutions  from  one  subset  is 
dominated  by  any  solution  in  the  subset.  These  conditions  ensure 
that  separable  operators  are  unambiguously  defined  in  the  sense  that 
any  two  distinct  optimally  ranked  bases  have  similar  structure 
(Lemma  2.2.1) . 

The  disjoint  subsets  of  an  optimally  ranked  basis  are  arranged 
in  linear  order  by  increasing  rate  of  growth.  All  possible  solutions 
of  an  inhomogeneous  difference  equation  with  arbitrary  right  side  are 
classifiable  by  comparison  with  the  solutions  in  an  optimally  ranked 
basis  of  the  corresponding  homogeneous  form.  Thus  a particular 
solution  y(r)  is  said  to  be  of  type  s if  (i)  y(r)  is  not  dominated  as 
r -r  CO  by  any  solution  in  the  first  s subsets  of  an  optimally  ranked 
basis;  (ii)  y(r)  is  dominated  by  every  solution  in  the  remaining 
a-s  subsets.  The  extreme  cases  s=a  and  s=0  correspond  to  particular 
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solutions  for  which  forward  recurrence  and  backward  recurrence  of  the 
given  difference  equation  are  stable,  respectively.  For  the  inter- 
mediate cases  0 < s < cr,  let  j be  the  number  of  solutions  in  the  first 
s subsets.  If  in  addition  to  being  separable  the  difference  operator  is 
(i, j ) -factorizable  for  some  initial  point  i,  analysis  of  the  boundary 
value  problem  shows  that  the  equivalent  pair  of  forward  and  backward 
recurrence  problems  (Theorem  1.2.2)  are  stable  in  the  sense  that  rounding 
errors  are  not  propagated  more  rapidly  than  the  solutions,  at  least  in 
the  asymptotic  region  belonging  to  the  given  problem.  Accordingly,  the 
solution  y(r)  of  the  boundary  value  problem  is  obtained  in  a stable 
manner  by  forward  elimination  without  pivoting  followed  by  back  sub- 
stitution (§3.3  and  Theorem  3.3.1).  Furthermore,  under  appropriate 
conditions  it  follows  that  (i)  y(r)  is  uniquely  determined  by  j initial 
values  (Theorem  3.1.1);  (ii)  the  algorithm  converges  pointwise  as  the 
terminal  point  proceeds  to  infinity  (Theorem  3.2.1). 

The  numerical  examples  included  in  this  thesis  were  chosen  so 
that  known  asymptotic  results  could  be  used  in  conjunction  with  the 
new  classification  theory  for  separable  operators.  This  enabled  the 
correct  number  of  initial  conditions  for  any  chosen  solution  to  be 
ascertained  immediately.  The  examples  also  show  that  it  is  sometimes 
desirable  to  associate  a solution  with  solutions  of  the  next  higher 
theoretical  type.  This  artificial  raising  of  the  type  is  manifested  in 
the  algorithm  by  the  specification  of  the  number  of  initial  conditions 
appropriate  for  the  higher  type.  This  artifice  preserves  unique  determination 


of  the  solution  and  also  convergence  of  the  algorithm,  but  introduces 
instability.  It  is  employed  only  when  the  asymptotic  separation  between 
two  t>n3es  is  weak,  consequently  the  loss  of  precision  resulting  from  the 
instability  is  mild.  The  compensating  advantage  is  that  the  rate  of  con- 
vergence is  increased  very  substantially,  enabling  a much  smaller  terminal 
point  to  be  used. 

The  trade-off  between  numerical  instability  and  slow  convergence  can 
be  regarded  heuristically  in  the  following  way.  First,  suppose  u(r)  ~ e 
and  v(r)  e are  two  fundamental  solutions  of  a second-order  equation. 

Then  one  initial  value  of  u(r)  determines  the  function  uniquely  and  the 
algorithm  is  both  convergent  and  stable.  Furthermore,  the  rate  of  convergence 
is  very  rapid.  This  can  be  seen  from  the  equation 


-n 

u (r)  = u(r)  - - — • v(r){l+o(l)}  , n , 

n n 

e 


where  solution  of  the  approximating  boundary  value  problem 

having  n as  its  terminal  point  (compare  Eq . 3.2.6).  The  rapid  convergence 
of  as  n ^ for  fixed  r is  readily  apparent.  Thus  the 

algorithm  performs  very  well  with  regard  to  both  numerical  stability  and 
rapidity  of  convergence  when  the  separation  of  solutions  is  strong  as  in 

this  example  of  exponential  separation. 

-1  -1/2 

Next,  suppose  that  u(r)  r and  v(r)  ~ r . Here  again,  one 

initial  value  of  u(r)  determines  the  function  uniquely  and  the  algorithm 
is  both  convergent  and  stable.  But  in  this  case 


u (r ) 
n 


-1 


u(r)  - 


-1/2 


• v(r){l+o(l)}  , 


n -f  00 
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and  it  is  apparent  that  the  rate  of  convergence  will  be  governed  by  the  very 

-1/2 

slowly  decreasing  function  n . Suppose,  however,  that  we  give  an 

additional  initial  value  of  u(r)  . With  two  initial  values  our  algorithm 

is  equivalent  to  forward  recurrence  of  the  given  difference  equation  and 

all  questions  of  convergence  disappear.  This  forward  recurrence  is,  of 

course,  unstable  for  the  desired  solution  u(r)  but  the  instability  is 

very  weak;  any  small  component  of  the  more  dominant  solution  v(r)  can 

1/2 

grow  no  faster  than  r relative  to  u(r)  . Although  this  rate  of  growth 

(or  any  other  positive  rate  of  growth)  is  ultimately  sufficient  to  obliterate 
the  desired  solution  when  computed  in  finite  precision  arithmetic,  subse- 
quences of  the  solution  of  considerable  length  can  be  computed  before  this 

happens.  Thus  in  many  practical  situations  where  the  separation  of  solutions 

T 

IS  weak,  a minor  modification  of  the  algorithm  performs  acceptably  well 

with  regard  to  both  stability  and  convergence. 

Finally,  suppose  that  the  relative  separation  of  u(r)  and  v(r)  is 
-3/2  -5/2 

like  r or  r (compare  Example  4.4.1  which  presents  the  computation 


of  a function  exhibiting  both  forms  of  behavior,  r 


-3/2 


on  the  subsequence 


-5/2 

of  odd  terms  and  r on  the  subsequence  of  even  terms).  These  separa- 

tion ratios  are  too  weak  to  allow  rapid  convergence  of  the  unmodified 
algorithm,  and  at  the  same  time  too  strong  for  the  numerical  instability  in 
the  modified  algorithm  to  be  truly  insignificant.  Computation  of  cases  such 
as  this  are  the  most  difficult  for  the  algorithm  to  perform  well,  and  per- 
formance can  be  improved  only  by  recourse  to  extended-precision  computation. 


In  general  the  modified  algorithm  will  require  a back  substitution  as  well 
as  a forward  elimination  stage.  Nevertheless,  similar  considerations  govern 
the  trade-off  between  numerical  instability  and  slow  convergence. 


Perhaps  it:  should  be  stressed  that  a complete  knowledge  of  the  asymp- 
totic behavior  of  all  of  the  solutions  of  the  homogeneous  equation  is 
essential  to  determine  the  correct  number  of  initial  conditions.  This  in- 
formation may  be  developed  for  a broad  class  of  linear  difference  equations 
using  known  analytical  results  (§2.2).  An  alternative  would  be  to  attempt 
to  construct  the  essential  classification  parameters  - the  number  of  dis- 
tinct types  of  solution,  and  the  number  of  linearly  independent  solutions 
of  each  type  - by  numerical  experimentation.  Under  this  approach,  the 
algorithm  presented  in  this  thesis  would  be  applied  with  chosen  initial  con- 
ditions appropriate  to  different  possible  realizations  of  the  classifica- 
tion theory,  that  is,  different  possible  forms  of  the  optimally  ranked  basis. 
Presumably,  incorrect  choices  would  be  revealed  by  discernible  numerical 
instability.  The  effectiveness  of  such  an  approach  has  not  be  assessed, 
however . 

In  conclusion,  the  results  presented  in  this  thesis  are  distinguished 
from  the  work  of  earlier  authors  in  the  following  ways.  First,  the  inter- 
plav  between  linear  difference  operators  and  infinite  band  matrices  achieves 
a fuller  blending  of  analytic  and  algebraic  ideas.  This  leads  to  a clearer 
understanding  of  stability,  including  the  tradeoff  choices  that  exist 
between  slow  convergence  and  mild  instability  when  the  asymptotic  separa- 
tion of  adjacent  solutions  is  weak.  Secondly,  the  classification  theory 
for  separable  operators  is  believed  to  be  new.  It  enables  the  proper 
number  of  initial  conditions  that  are  prescribed  by  the  convergence 
theorem  to  be  determined  directly  from  the  asymptotic  behavior  of  fundamen- 
tal solutions  of  the  given  difference  equation.  In  contrast,  Zahar  employs 
fundamental  solutions  of  the  adjoint  equation  and  Oliver  applies  3 
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point-by-point  growth  criterion  between  pairs  of  solutions.  Finally,  the 
numerical  examples  presented  here  have  been  chosen  in  such  a way  that  the 
solutions  exhibit  a much  richer  variety  of  growth  rates  than  in  any  pre- 


vious investigation. 


APPENDIX 


Proof  of  Lemma  2.1.1 

The  proof  consists  of  ( i)  identifying  the  sequence  defined  by 
equation  (2.1.4)  as  an  almost  periodic  sequence,  and  (ii)  proving  that 
evem,’  ncr-zero  almost  periodic  sequence  has  an  infinite  subsequence  which  is 
bounded  away  from  zero.  L'e  begin  by  presenting  the  essential  elements  of 
the  theoiq/  of  almost  periodic  functions. 

First  we  give  the  characterization  of  almost  periodic  functions  that 
is  taken  as  the  definition  by  Corduneanu  [7  ].  A function  of  the  form 

T(x)  = I c e , - » < X < “ , 

k=l 

where  the  c,  are  complex  numbers,  the  \ are  real  numbers  and 
.<  ■ k, 

i = V-1  , is  called  a complex  trigonometric  polynomial . A function  f (x) 
taking  complex  values  and  having  all  real  numbers  as  domain  is  an  almost 
periodic  function  provided  that  to  each  s > 0 there  corresponds  some 
trigonometric  polynomial  T_ (x)  such  that 

if(x)  - T^(x) 1 < e 

uniformly  on  - » < x < » . Clearly,  every  complex  trigonometric  poly- 
nomial is  itself  an  almost  periodic  function. 

An  alternative  characterization  is  the  definition  of  H.  Bohr  [4  ]: 

A continuous  complex  valued  function  f(x)  defined  on  the  real  line  is 
almost  periodic  if  to  each  z > 0 there  corresponds  a positive  real 
number  ^ such  that  every  open  interval  of  the  real  line  of  length 

contains  at  least  one  point  ^ for  which 


I f (x+O  - f (x)  I < e 
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uniformly  for  every  x . Any  number  ^ which  satisfies  this  condition 
is  called  a translation  number  of  f corresponding  to  s , or  more  brief- 
ly, an  5 -translation  number  of  ^ , 

An  t-translation  number  of  f is  "almost"  a period  of  f . Every 
periodic  function  is  almost  periodic  as  well.  If  p is  the  fundamental 
period  of  a periodic  function  f , then  rp  for  r = 0,+l,+2,...  are 
translation  numbers  of  f for  each  s . Bohr's  definition  generalizes 
the  fact  that  the  periods  of  a periodic  function  form  an  arithmetic 
progression. 

Bohr's  definition  and  Corduneanu's  definition  define  exactly  the 
same  class  of  functions.  This  is  the  "fundamental  theorem"  of  the  theory 
of  almost  periodic  functions.  See  either  [4  ] or  [7  ] for  a full  account. 

Turning  now  to  sequences,  we  have  the  following  definition  presented 
in  [ 7 ,§1.6]:  A function  f(t)  , taking  complex  values  and  having  the  set 

of  integers  as  domain,  is  an  almost  period  sequence  provided  that  to 
each  e > 0 there  corresponds  some  positive  integer  N = N(s)  such  that, 
within  any  set  of  N consecutive  integers,  there  exists  an 
integer  p for  which 

, f (r+p)  - f (r)  ! < £ . 

uniformly  over  all  integers  r . Corduneanu  proves  "A  necessary  and 
sufficient  condition  for  a sequence  almost  periodic  is  the 

existence  of  an  almost  periodic  function  f(x)  , - =»  < x < “=  , such  that 
= f(r)  , r = 0 , j;^l  ,+2,  . . . " ; see  [7,  Theorem  1.27]  . 

VTe  are  now  ready  to  prove  Lemma  2.1.1.  If  r were  allowed  to  assume 
all  real  values,  instead  of  being  restricted  to  r = 0,1,2,...,  then 
equation  (2.1.4)  defines  a trigcnometric  polynomial.  By  Corduneanu's 
definition,  every  trigonometric  polynomial  is  an  almost  periodic  function. 


By  the  theorem  quoted  above,  the  sequence  x(r)  , r = 0,+l,+2,... 
defined  by  equation  (2.1.4)  is  an  almost  periodic  sequence.  The  only 
thing  remaining  is  to  prove  the  following  assertion: 

If  an  almost  periodic  sequence  {f(r)}  is  not  identically  zero, 
then  there  exists  an  infinite  subsequence  of  {f(^)}  which  is  bounded 
away  from  zero. 

Proof  Assume  {f(r)}  is  not  identically  zero.  Let  be  such  that 

fCr^)  ^ 0 and  put  A = |f(rQ)|  . Put  e = A/(2+A)  . Then  e $ (0,1)  , 

and  by  the  definition  of  an  almost  periodic  sequence  there  exists  a 
positive  integer  , say,  such  that  for  some  integer  p^  c {1,2,...,N^} 

the  inequality 


f (^o+Pi)  - f (^q)  I < e 


is  satisfied.  Put  r^  = • Next,  there  exists  a positive  integer 

such  that  for  some  integer  p^  - {1,2,...,N2}  the  inequality 


f (ti+Pi)  - f (t^)  I < £ 


is  satisfied.  Put  ^^2  ~ ^1  ^2  ' Continuing  in  this  manner,  we  see  that 

for  each  k = 3,4,5,...  we  can  find  an  integer  ^ ]_  such  that 


1,2,3,...  we  have 


f(r^)  - f(r^j)|  S |f(r.)  - 


r j £ 1 . 

< 4 & < irr  < T ^ 

j=i 


For  each  k = 


Since 


A = ! f (i^g)  I . we 

!f(r. 


have 

^)|  = |f(rg)  + - f(rg)| 

2 A - |f(r^)  - f(rg)|  > 

• m 


which  finishes  the  proof  of  the  assertion 


N)  IH" 
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